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Abstract 


The gas foil bearing (GFB) technology is considered one of the key factors for the 
intended transition to oil-free rotating machinery in future transportation systems. 
Besides numerous advantages in terms of size, weight, efficiency, and cleanliness, 
GFBs offer the unique ability to be lubricated with working fluids such as refrigerants. 
However, the computational analysis of refrigerant-lubricated GFB-rotor systems 
represents an interdisciplinary problem of enormous complexity and with conflicting 
interests between all-encompassing but efficient modeling and solution approaches. 
This thesis succeeds in exploring and pushing forward existing limits of feasibility 
and thereby establishes a new strategy that enables stability and bifurcation analyses. 


Owing to the precisely identifiable fluid-structure-rotor interaction mechanisms, 
three submodels that are of reasonable complexity but nevertheless take into account 
all relevant nonlinearities can finally be transformed into a single dynamical system. 
As an essential model feature, the non-ideal characteristics of a typical refrigerant 
that may undergo vapor-liquid phase transitions are described by thermodynamic 
equations of state to be included in modified Reynolds and temperature equations. 
Also, it becomes feasible with the proposed lumped-element foil structure model 
to include dry friction in various ways reaching from highly efficient regularized 
Coulomb friction models to bristle models that can reproduce stick-slip transitions. 
Altogether, this makes the entire problem accessible to rigorous mathematical theory 
and allows for developing a monolithic research code with interchangeable modules. 


Introduced by a brief overview of the mathematical theory of dynamical systems, 
the description and discussion of relevant findings is divided into three sections and 
guides through a selection of the most important nonlinear effects and phenomena. 
Firstly, the particular features of refrigerant lubrication with phase transitions are 
focused on with regard to both steady-state operation and vibrational self-excitation. 
Secondly, the influence of some foil structure design parameters and the importance 
of dry friction are thoroughly investigated using stability and bifurcation analyses. 
Finally, the addition of rotor unbalance gives an outlook to quasi-periodic behavior 
with even more complex scenarios resulting from combined excitation mechanisms. 


Keywords: Aeronautics, Bifurcation Analysis, Dry Friction, Fluid Film Lubrication, 
Gas Foil Bearing, Rotor Dynamics, Thermodynamics, Tribology, Vapor Cycle System. 
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Erweiterte Kurzfassung 


Die nachfolgende Kurzfassung in deutscher Sprache bietet einen ersten Uberblick 
uber die wesentlichen Methoden sowie Erkenntnisse der vorliegenden Dissertation. 


Einleitung und Motivation 


Bis zum heutigen Tag stellt das Auftreten schädlicher oder sogar zerstörerischer 
Vibrationen ein ernstzunehmendes Problem in all jenen Ingenieurdisziplinen dar, 
welche sich mit rotierenden Maschinen befassen. Neben zumeist unvermeidlichen 
drehzahlsynchronen Schwingungen aufgrund der Restunwucht von Rotoren können 
Fluidlager darüber hinaus selbsterregte subsynchrone Schwingungen verursachen. 
Trotz dieses bedeutsamen Nachteils, der schon aus den Anfängen der Ölschmierung 
bekannt ist, entwickelt sich neuerdings reges wissenschaftliches Interesse vor allem 
an Gasfolienlagern (GFBs'). Indem sie auf demselben Prinzip der fluiddynamischen 
Schmierung beruhen wie einfache Gleitlager, funktionieren solche GFBs selbsttätig 
und verlassen sich auf die Ausbildung eines Schmierkeils, wodurch viskoses Fluid 
in einen sich verengenden Spalt hineingezogen und somit Druck aufgebaut wird. 
Sobald die Welle eines schnelldrehenden Rotors vollständig von der umgebenden 
Struktur getrennt und in ein durchgängiges Fluidpolster eingebettet ist, führt das 
Nichtvorhandensein von Festkörper-zu-Festkörper-Kontakt zu erstaunlich kleinen 
Lagermomenten und ermöglicht dabei einen weitestgehend verschleißfreien Betrieb. 


Im Vergleich zu flüssigkeitsgeschmierten Lagern reduziert die um ein Vielfaches 
niedrigere Viskosität von Gasen die Verlustleistung noch weiter, erfordert jedoch auch 
deutlich höhere Drehzahlen bis zum Aufbau eines lasttragenden Schmierfilms und 
macht das System anfälliger für Instabilitäten. Als Abhilfemaßnahme hinsichtlich der 
nur schwach ausgeprägten Fluiddämpfung bringt eine nachgiebige Folienstruktur 
trockene Reibung in das System ein und schwächt so selbsterregte Schwingungen ab. 
Insgesamt ergibt dies eine Technologie mit außergewöhnlicher Leistungsfähigkeit 


t Das verwendete Akronym GFB ist abgeleitet von der englischen Bezeichnung Gas Foil Bearing. 
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Erweiterte Kurzfassung 
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(a) Technische Zeichnung einer AiResearch-VCM aus den 
1980ern mit kaltemittelgeschmierten GFBs, zuerst verbaut in 
Lockheed P-3C der US Navy (in Anlehnung an [DWHg0)). 


(b) Fotografie eines Open-Source-GFBs 
(Bump-Type-GFB der ersten Generation) 
von der NASA (Ausschnitt aus [NASo7]). 


Abbildung 1: Anschauliche Beispiele für GFBs in Anwendungen aus der Luftfahrt. 


und Zuverlässigkeit, selbst unter so anspruchsvollen Betriebsbedingungen wie etwa 
extremen Temperaturverhältnissen. Für eine Vielzahl von Anwendungen erwächst 
der herausragendste Vorteil von GFBs jedoch aus der Einsetzbarkeit grundsätzlich 
beliebiger Arbeitsmedien als Schmiermittel, was den Weg zu ölfreien rotierenden 
Maschinen bereitet, die jegliches Risiko von Umgebungskontamination ausschließen. 


Die Entwicklung von GFBs steht in enger Verbindung zu deren häufiger Verwendung 
in Klimaanlagen von Flugzeugen, bei welchen zwischen Luftkältemaschinen (ACMs?) 
und wie in Abb. 1a gezeigten Dampfkältemaschinen (VCMSs5) zu unterscheiden ist. 
In VCMs kommen statt der in ACMs genutzten Umgebungsluft spezielle Kältemittel 
wie zum Beispiel der Halogenkohlenwasserstoff R-245fa zum Einsatz, um im Rahmen 
eines thermodynamischen Kreisprozesses ständige Phasenübergänge zu durchlaufen. 
Als mögliche Umsetzung mit zweiphasigem R-245fa schmierbarer Lager stellt die 
Fotografie in Abb. 1b die weitestverbreitete Variante eines Bump-Type-GFBs der 
ersten Generation dar, welche oft auch mit mehreren Folienpads konstruiert werden. 


Die Untersuchung in GFBs gelagerter Rotoren setzt Kenntnisse in verschiedensten 
Bereichen voraus, darunter etwa Rotordynamik, Fluiddynamik, Thermodynamik, 
Strukturdynamik oder Tribologie. Dementsprechend stellt sich die weltweit laufende 
Forschung als interdisziplinäre Aufgabe dar, welche gegenwärtig immer stärker auf 
rechnergestützte Werkzeuge setzt. Nach wie vor bedeutet die schiere Komplexität 


? Das verwendete Akronym ACM ist abgeleitet von der englischen Bezeichnung Air Cycle Machine. 
3 Das verwendete Akronym VCM ist abgeleitet von der englischen Bezeichnung Vapor Cycle Machine. 


Modellierung von GFB-Rotor-Systemen 


des Problems jedoch eine erhebliche Erschwernis für das tatsächliche Erreichen eines 
allumfassenden und dennoch effizienten Modellierungs- wie auch Lösungsansatzes. 
Insbesondere mangelt es der Literatur an Minimalmodellen kältemittelgeschmierter 
GFB-Rotor-Systeme, die einer monolithischen Herangehensweise zugänglich wären. 
Angesichts dieser Forschungslücke setzt sich die vorliegende Dissertation das Ziel 
bisherige Machbarkeitsgrenzen zu überwinden, sodass erstmals rotordynamische 
Analysen zum Systemverständnis mit GFBs ausgestatteter VCMs beitragen können. 


Modellierung von GFB-Rotor-Systemen 


Aufbauend auf die Festlegung der Lagergeometrie und Ermittlung grundlegender 
kinematischer Zusammenhänge schafft Kapitel 2 nötige Voraussetzungen für die 
letztliche Aufstellung eines gekoppelten Fluid-Struktur-Rotor-Interaktionsproblems. 


Die Skizze in Abb. 2a zeigt einen Querschnitt entlang der Mittellinie eines GFBs mit 
drei Folienpads sowie den darin montierten Wellenzapfen eines zugehörigen Rotors. 
Zur Veranschaulichung der relevanten kinematischen Verhältnisse übertreibt die 
noch weiter abstrahierte Darstellung in Abb. 2b die Abmessungen des Lagerspiels C 
deutlich gegenüber dem Radius R. Als weitere wichtige geometrische Parameter 
werden die Anzahl Np und auch die Ausrichtung xı der Pads eingeführt, sodass die 
gesamte Modellierung möglichst universell alle denkbaren Konfigurationen abdeckt. 


Zur Beschreibung von Feldgrößen im Schmierspalt werden die Koordinaten p und z 
benötigt, wobei stets auch die Abhängigkeit von der Zeit t Berücksichtigung findet. 
In möglichst allgemeiner Darstellung ist das Verformungsfeld der Folienstruktur 
dann durch q(¢,t) gegeben, wohingegen die Verlagerung eines Zapfens durch £(t) 
und y(t) oder ebenso durch die Exzentrizität e(t) und die Winkellage T (t) feststeht. 
Die Rotordrehzahl no und damit die Kreisfrequenz wg = 27n0 werden grundsätzlich 
als konstant angesehen, sodass Hoch- und Auslaufvorgänge quasistationär erfolgen. 


Aus einer trigonometrischen Betrachtung am hervorgehobenen Dreieck in Abb. 2b 
folgen allgemeingültige Ausdrücke für die in der Exzentrizität linearisierte Filmdicke 


h(g,t A — e(t) cos[p — T(#)] — q(t) (1a) 
— [&(#) sin g + y(t) cos 9] — q(t) (1b) 

T t) — q(p,t) (Go) 

als Uberlagerung des Lagerspiels, des sinusférmigen Beitrags der Zapfenverlagerung 


und des Verformungsfelds der Folienstruktur. Zur vollstandigen Kopplung zwischen 
Fluid, Struktur und Rotor wird zudem die Zeitableitung von Gleichung (1) gebildet. 
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Erweiterte Kurzfassung 


(a) Darstellung inspiriert vom typischen Design der (b) Abstrakte Geometrie mit übertriebenem Spalt 
gängigen Bump-Type-GFBs der ersten Generation. zur Veranschaulichung kinematischer Verhältnisse. 


Abbildung 2: Querschnittzeichnungen eines beispielhaften Drei-Pad-GFBs mit eingesetzter Rotorwelle, 
die zeigen wie das Fluid (in Blau) den rotierenden Zapfen (in Rot) von der umliegenden Struktur trennt. 
Die Pads bestehen aus glatten Deckfolien (in Dunkelgrün) gestützt durch gewellte Tragfolien (in Hellgrün), 
welche sich unter dem Einfluss des strömungsinduzierten Drucks verformen und die Spaltform festlegen. 


Letztendlich kann ein modifiziertes Jeffcott-Laval-Rotormodell aufgestellt werden, 
in welchem die Masse nicht vollständig in der Scheibe konzentriert, sondern zum 
Teil in die Zapfen verschoben ist. Unter Berücksichtigung des Eigengewichts sowie 
einer statischen Unwucht folgt so ein System gewöhnlicher Differentialgleichungen, 
worin die Lagerkräfte der GFBs als Integrale über ein Druckfeld p(,z,t) eingehen. 
Mit der beschriebenen Herangehensweise wäre auch eine Betrachtung komplexerer 
Rotormodelle möglich, wie sie etwa für praktische Auslegungen oft benötigt werden. 


Schmierfilmtheorie für phasenveränderliche Fluide 


Die bekannte Schmierfilmtheorie erlaubt eine Berechnung von Druckverteilungen in 
engen Spalten auf Grundlage vereinfachter Gleichungen der Stromungsmechanik. 
Eine Betrachtung kältemittelgeschmierter GFBs stellt jedoch einige der üblicherweise 
getroffenen Annahmen infrage, weshalb Kapitel 3 die Theorie sorgfältig aufarbeitet. 


Wendet man die Gesetze der Gleichgewichtsthermodynamik auf ein Kontinuum an, 
ist jeder mögliche Zustand des Schmiermittels durch Vorgabe der Massendichte p 
und der Temperatur # lokal eindeutig beschrieben. Für ein gegebenes Kältemittel 
lassen sich damit insbesondere alle weiteren thermodynamischen Zustandsgrößen 
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Schmierfilmtheorie fiir phasenveränderliche Fluide 


als davon abhängige Zustandsgleichungen formulieren, darunter etwa der Druck p, 
die spezifische Warmekapazitat cy, der Kompressionskoeffizient ay oder auch die 
dynamische Viskosität u. Im interessierenden Betriebsbereich zeigen sich für R-245fa 
besonders ausgeprägte Nichtlinearitäten beim Überschreiten einer Phasengrenzlinie, 
wodurch sich das ideale Gasgesetz etwa im Zweiphasengebiet als ungeeignet erweist. 


Die Geschwindigkeitsprofile von Strömungen in dünnen Schmierfilmen lassen sich 
durch vereinfachte Navier-Stokes-Gleichungen gut annähern und führen über die 
Kontinuitätsgleichung zur modifizierten Reynolds-Gleichung für nicht-ideale Gase 


f] ‚wo ð _ 1 à [ ph? ap ‚9 ph? op 
ot (er) ' 2 ap (er) R209 (= op) ` əz \ 12u dz J` @) 
Diese partielle Differentialgleichung beschreibt direkt den Druckaufbau im Fluid, 


muss aber zunächst mittels thermodynamischer Zustandsgleichungen geschlossen 
und durch das Einsetzen von Gleichung (1) an Struktur und Rotor gekoppelt werden. 


Darüber hinaus erweist es sich wegen der temperaturempfindlichen Eigenschaften 
von Kältemitteln als unabdingbar, die Energiegleichung für den dünnen Schmierfilm 
gleichermaßen zu vereinfachen zu einer Temperaturgleichung für nicht-ideale Gase 


op | 98 , 2098 _ h? ( 1 ap ad | dpae 
| Te 2 ap 12u \ R20p0@ ` dz az 


gpn wh 1 2 (I ap) (ap 
OP oe 2 og R2d9\ 12u ðP dz | 12y dz 
u h | 1 [op i op | 
L R2,.52 | | 
Bo( 8 - ds) + Ru? aor (3) ! (2) 6) 


Wie in die Reynolds-Gleichung fließen hier Zustandsgleichungen und Filmdicke ein, 
sowie zudem der Wärmeübergangskoeffizient fọ und die Umgebungstemperatur Oo. 


Als stark nichtlineare partielle Differentialgleichungen besitzen sowohl Gleichung (2) 
als auch Gleichung (3) gemischt parabolisch-hyperbolischen Charakter, was einer 
Wahl physikalisch sinnvoller Randbedingungen gewisse Anforderungen auferlegt. 
Nachdem die Gebiete, auf denen nach Lösungen gesucht wird, samt zugehöriger 
Ränder identifiziert sind, müssen einige Annahmen über den Zustand des die GFBs 
umgebenden Fluids getroffen werden, woraus Dirichlet-Randbedingungen folgen. 
Insgesamt ergeben sich so als Fluidmodell mehrere zu lösende Randwertprobleme, 
deren Anzahl von einer den Folienpads entsprechenden Gebietsaufteilung abhängt. 
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Erweiterte Kurzfassung 


fia Y 


— 


FOR RORNOR Ya Z 
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X1,1 f X1,2 q x13 f X14 q X15 
j<— R 


Abbildung 3: Skizze des Modells mit konzentrierten Parametern ftir ein Pad mit fiinf Stiitzelementen, 
das die Tragfolie („bump foil” in Hellgrün) auf eine Anordnung von Stäben, Federn und Massen reduziert, 
während die Deckfolie („top foil” in Dunkelgrün) linear zwischen den Stützelementen interpoliert wird. 


Effiziente Folienmodelle mit Reibungseigenschaften 


Da die Folienstruktur von entscheidender Bedeutung für das Systemverhalten ist, 
existieren in der Literatur bereits eine Vielzahl entsprechender Modellierungsansätze, 
wovon die meisten jedoch übertrieben detailliert und dadurch recht ineffizient sind. 


Die Skizze in Abb. 3 zeigt das neue Minimalmodell mit konzentrierten Parametern 
aus Kapitel 4, das den Aufbau eines Tragfolienabschnitts nachempfindet, indem jedes 
der Stützelemente einer Anordnung zweier starrer Stäbe und einer Feder entspricht. 
Ebenso wie die Steifigkeiten kg werden auch die den Stützelementen zugeordneten 
Punktmassen mp anhand der Eigenschaften industriell eingesetzter Folien bestimmt. 
Auch wenn Eigenfrequenzen der Folien für die Systemdynamik irrelevant scheinen, 
erweist sich eine Mitnahme der Trägheit als zweckmäßig für numerische Methoden. 


Zur kinematischen Beschreibung des Modells werden generalisierte Koordinaten Xn,m 
als Verschiebungen an den Punktmassen eingeführt, was im Hinblick auf die dort zu 
berücksichtigende Reibung besonders günstig ist. Schließlich gelingt es außerdem, 
auch die Verschiebungen yy, der Spitzen auf ebendiese Darstellung zurückzuführen, 
wofür sich auf Geschwindigkeitsebene eine Jacobi-Matrix als zielführend erweist. 
Zwischen den Stützelementen wird die Deckfolie einfach durch lineare Interpolation 
angenähert und liefert damit das Verformungsfeld als Begrenzung des Schmierspalts. 


Durch die regelmäßige Geometrie des Systems und die vielen Zwangsreaktionen 
können die Bewegungsgleichungen über Lagrangesche Gleichungen zweiter Art 
mühelos hergeleitet werden. Es ergibt sich so für jedes Pad ein System gewöhnlicher 
Differentialgleichungen von gleicher Form wie für eine harmonische Schwingerkette. 
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Numerische Analyse von GFB-Rotor-Systemen 


Die von oben wirkenden Kräfte F,,, entsprechen der Belastung durch das Fluid und 
resultieren aus der Integration über das Druckfeld in einer zugehörigen Umgebung. 


Die Reibkräfte mit Koeffizient ug in Kontakten zwischen Tragfolie und Lagerschale 
sind für Haften wie auch für Gleiten gegeben durch den mengenwertigen Ausdruck 


{—ABFn,m } falls Xnm < 0, 
fam € [-HBFum, +ugBFn,m] falls Žn,m =0, (4) 
{+upEnm } falls Xnm > 0. 


Da Reibung in Form von Gleichung (4) numerisch sehr schwer zu beherrschen ist, 
kommt man kaum umhin, ausgeprägtes Haften zu vernachlässigen und etwa mittels 


x 
fam = UBEn,m tanh (ur) (5) 


regularisierte Coulombsche Reibung einzuführen oder alternativ auf ein geeignetes 
Bürstenmodell zurückzugreifen, welches Haft-Gleit-Übergänge abzubilden vermag. 


Numerische Analyse von GFB-Rotor-Systemen 


Bevor das Gesamtsystem rechnerisch untersucht wird, erscheint eine einheitliche 
Entdimensionierung aller aus der Modellierung erhaltener Gleichungen vorteilhaft. 
Im Zuge dieses Vorgehens tauchen eine ganze Reihe physikalisch interpretierbarer 
dimensionsloser Kennzahlen auf, welche jedes Problem anschaulich charakterisieren. 


Zur näherungsweisen Lösung der Reynolds- und Temperaturgleichung werden die 
räumlichen Ableitungen anhand eines Finite-Differenzen-Verfahrens diskretisiert, 
wobei sich ein Upwind-Schema als am besten geeignet erweist. Schlussendlich 
überführt dies das Fluidmodell in ein System gewöhnlicher Differentialgleichungen. 


Die in dieser Arbeit vorgeschlagene Modellierung erlaubt es, den Systemzustand 

jederzeit eindeutig zu beschreiben durch einen Fluid-Struktur-Rotor-Zustandsvektor 
T 

s=| si st sk |enck, (6) 

der zu einem hochdimensionalen Phasenraum (2 gehört. Nach kurzer Umformung 

aller zuvor hergeleiteten Modellgleichungen ergibt sich so das dynamische System 


D = R{S, Ao}, R:QxR->RN, (7) 
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Erweiterte Kurzfassung 


Yeo 
XR ~! i: 7 \ 
| 
e=0 : |. a} 
Ag = 0,2\ |” 
Net 
(a) Fixpunktlösung. (b) Periodische Lésung. (c) Quasiperiodische Lésung. 


Abbildung 4: Plots transienter Bahnkurven der Scheibe des Rotors aus numerischer Zeitintegration, 
welche einen ersten Eindruck vom dynamischen System langfristig erreichbarer Lösungstypen verschaffen. 


worin die dimensionslose Drehzahl oder Lagerzahl A5 Bifurkationsparameter ist. 
Obwohl also die Teilmodelle für Fluid, Struktur und Rotor separat entworfen sind, 
führt dieser entscheidende Schritt zu einer streng monolithischen Lösungsstrategie. 
Basierend auf Kapitel 5 wird der Forschungscode GFB++ in C/C++ entwickelt, 
welcher von einer hohen Ausführungsgeschwindigkeit profitiert und gleichzeitig 
dank objektorientierten Programmierparadigmas modular konzipiert werden kann. 


Nichtlineare Effekte in GFB-Rotor-Systemen 


In Abhängigkeit der vielen Parameter kann sich das nichtlineare dynamische System 
auf lange Sicht sehr unterschiedlich verhalten, woraus sich gemäß den Plots in Abb. 4 
eine Zuordnung zu verschiedenen Lösungstypen ergibt. Erschwerend kommt hinzu, 
dass Lösungen instabil werden können und oft mit weiteren Lösungen koexistieren, 
was mit Methoden wie etwa einer Brute-Force-Zeitintegration kaum zu erfassen ist. 


Eine bewährte Methode zur systematischen Untersuchung dynamischer Systeme 
ist die Lösungsverfolgung mit gleichzeitiger Stabilitäts- und Bifurkationsanalyse. 
Stellvertretend für die lange Reihe damit durchgeführter Studien zeigen die Plots in 
Abb. 5 den Einfluss des dimensionslosen Rotordämpfungsparameters Dr auf das von 
der Lagerzahl Ag abhängige Verhalten kältemittelgeschmierter GFB-Rotor-Systeme. 
Offensichtlich unterscheiden sich die Bifurkationsszenarien nach Abb. 5a deutlich, 
sodass entweder eine unter- oder eine überkritische Hopf-Bifurkation (HB) auftritt, 
der wiederum bis zu zwei Sattelpunkte von Grenzzyklen (LPCs*) folgen können. 


4 Das verwendete Akronym LPC ist abgeleitet von der englischen Bezeichnung Limit Point of Cycles. 
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Nichtlineare Effekte in GFB-Rotor-Systemen 


Einhüllende min Yk- max Yr 


Verhältnis Asub / Ao 


15 T T T T T Dämpfung Dr 
10 + —— 0,00 x 1072 
i e —— 150x10 
ag weg = | «| — 3235x107 
0,0 
—0,5 
—1,0 
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Dimensionslose Drehzahl Ao 


(a) Bifurkationsdiagramme für Rotorscheibenauslenkungen (oben) und Frequenzverhältnisse (unten). 


(b) Ausgewählte Bahnkurven der Rotorscheibe (obere Reihe) und der Wellenzapfen (untere Reihe). 


Abbildung 5: Ergebnisse der Lösungsverfolgung und Bifurkationsanalyse für unterschiedliche Dämpfung 
(durchgezogen: stabil, gestrichelt: instabil, Markersymbole: Bifurkationen bzw. verzweigende Lösungen). 
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Erweiterte Kurzfassung 
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A a Umfangskoordinate p ne Umfangskoordinate p 
(a) Dampfanteilverteilung (Op = 144, Np = 3). (b) Dampfanteilverteilung (©, = 144, Np = 1). 


Abbildung 6: Plots zur Veranschaulichung von Kondensation in Drei-Pad- und Ein-Pad-GFBs. 


Auf praktische Anwendungen übertragen, verursacht das Überschreiten einer HB 
den Stabilitätsverlust eines stationären Betriebspunkts, woraufhin sich unerwünschte 
selbsterregte Schwingungen einstellen. Während sich aus dem oberen Plot in Abb. 5a 
die Amplituden ablesen lassen, liefert der untere Plot die zugehörigen Frequenzen 
bezogen auf Ag, sodass der Wert 1⁄2 dem typischen Halbfrequenzwirbel entspräche. 
In Abb. 5b sind schließlich sowohl statische Verlagerungsbahnen als auch Orbits von 
Grenzzyklen projiziert auf Rotorscheiben- und Wellenzapfenkoordinaten dargestellt. 


Eine Besonderheit aller kältemittelgeschmierten GFBs besteht darin, dass im Betrieb 
grundsätzlich Phasenübergänge zwischen gasförmig und flüssig auftreten können, 
was in Abhängigkeit von Parameterwerten bei allen Lösungstypen beobachtet wird. 
Die Plots in Abb. 6 zeigen den berechneten Dampfanteil als über den Schmierspalt 
verteilte Feldgröße in Drei-Pad- und Ein-Pad-GFBs an stationären Betriebspunkten, 
wobei der dimensionslose Wärmemanagementparameter ©, = 144 so gewählt ist, 
dass sich das Fluid nicht übermäßig erhitzt. Infolge lokaler Kondensation verändert 
sich beispielsweise auch die zugehörige Druckverteilung deutlich in ihrer Gestalt, 
was in allerletzter Konsequenz der Tragfähigkeit derartiger GFBs zum Nachteil gerät. 


Neben Betrachtungen zur Rolle trockener Reibung und zum Einfluss einer Unwucht, 
werden im weiteren Verlauf von Kapitel 6 außerdem Parameterstudien durchgeführt, 
deren zentrale Erkenntnisse sich in nachfolgender Zusammenfassung widerspiegeln. 


Zusammenfassung und Ausblick 


Die vorliegende Dissertation befasst sich thematisch mit der Dynamik von Rotoren 
in kältemittelgeschmierten GFBs und erarbeitet einen neuen Modellierungsansatz 
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Zusammenfassung und Ausblick 


zur Ermöglichung von Stabilitäts- und Bifurkationsanalysen. Dank der präzise 
erfassbaren Fluid-Struktur-Rotor-Interaktionsmechanismen können dabei letztlich 
drei Teilmodelle, welche trotz Berücksichtigung aller relevanten Nichtlinearitäten 
recht überschaubar bleiben, in ein einziges dynamisches System überführt werden. 
Als wesentliches Modellmerkmal beschreiben in Reynolds- und Temperaturgleichung 
einbezogene thermodynamische Zustandsgleichungen das nicht-ideale Verhalten 
eines typischen Kältemittels mit Phasenübergängen zwischen gasförmig und flüssig. 
Darüber hinaus erlaubt das vorgeschlagene Folienstrukturmodell mit konzentrierten 
Parametern die Einbindung trockener Reibung auf mehrere verschiedene Weisen, 
von effizienter regularisierter Coulombscher Reibung bis hin zu Bürstenmodellen mit 
Erfassung von Haft-Gleit-Übergängen. Insgesamt macht dies das Gesamtproblem der 
strengen mathematischen Theorie dynamischer Systeme zugänglich und ermöglicht 
den Aufbau eines monolithischen Forschungscodes mit auswechselbaren Modulen. 


Bei Betrachtung des Langzeitverhaltens kältemittelgeschmierter GFB-Rotor-Systeme 
entwickeln sich transiente Zustände in Abhängigkeit von Auslegungsparametern 
sowie der systematisch variierten Rotordrehzahl zu Fixpunktlösungen, periodischen, 
quasiperiodischen oder chaotischen Lösungen. Für perfekt ausgewuchtete Rotoren 
wird dabei beobachtet, dass Gleichgewichtslagen ihre Stabilität bei entsprechender 
kritischer Drehzahl mittels unter- oder überkritischer Hopf-Bifurkation verlieren, 
wodurch ganz typischerweise unerwünschte selbsterregte Schwingungen entstehen. 
Unter dann vorliegender dynamischer Belastung zeigen sich kältemittelgeschmierte 
GFBs besonders anfällig für Phasenübergänge, was schließlich zum Zusammenbruch 
des Schmierfilms und zum Systemversagen führen kann. Im stationären Betrieb wird 
Kondensation hingegen eher bei niedrigen Drehzahlen zum Problem, obwohl sich 
GFBs während des Betriebs üblicherweise bis zu einem solchen Ausmaß erhitzen, 
dass der Sättigungsdampfdruck nirgendwo erreicht wird. Um GFBs im Allgemeinen 
und oft zu Lasten der Tragfähigkeit weniger schwingungsempfindlich zu machen, 
entpuppen sich Anzahl und Ausrichtung der Pads als einflussreichste Eigenschaften. 
Außerdem lässt sich zeigen, dass die Nachgiebigkeit und Beweglichkeit einzelner 
Folienabschnitte in Verbindung mit der Reibung eine entscheidende Rolle spielen, 
um die Schwingungsamplituden überhaupt auf ein verträgliches Maß zu reduzieren. 
Am Ende führt das Hinzunehmen einer realistischen Rotorunwucht keineswegs zu 
einer Entkräftung vorheriger Erkenntnisse, sondern vielmehr zu deren Erweiterung 
hin zu noch komplexeren Szenarien aufgrund kombinierter Anregungsmechanismen. 


Schlussendlich werfen die Erkenntnisse der vorliegenden Dissertation eine Reihe von 
Folgefragen auf und liefern zahlreiche Anknüpfungspunkte für zukünftige Arbeiten, 
beispielsweise im Hinblick auf Modellverfeinerungen oder experimentelle Versuche. 


Die Dissertation schließt mit einem Symbolverzeichnis sowie einer Literaturübersicht. 
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Résumé étendu 


Le résumé suivant en langue française constitue un premier aperçu des principales 
méthodes utilisées ainsi qu’un bilan des résultats essentiels de ce mémoire de these. 


Introduction et motivation 


Jusqu’a ce jour, l'apparition de vibrations néfastes ou même destructrices représente 
un problème sérieux dans toutes les disciplines d'ingénierie qui s'intéressent aux 
machines tournantes. Outre les oscillations synchrones avec la rotation, qui sont 
généralement inévitables en raison du balourd résiduel d’un rotor, les paliers à film 
fluide peuvent également induire des vibrations sous-synchrones et auto-entretenues. 
Malgré cet inconvénient important, déjà connu depuis les débuts de la lubrification 
à l'huile, un énorme intérêt scientifique se manifeste récemment en particulier 
pour les paliers à feuilles (GFBs*). Sur la base du même principe de lubrification 
dynamique que les simples paliers lisses, de tels GFBs opèrent automatiquement et 
reposent sur la formation d’un coin de fluide dans lequel un lubrifiant visqueux est 
entraîné dans un espace convergent, ce qui finalement fait augmenter la pression. 
Dès que l'arbre d’un rotor à grande vitesse est entièrement séparé de la structure 
environnante et qu'il flotte sur un coussin de fluide continu, l'absence de contact 
solide sur solide se traduit par des couples de frottement particulièrement faibles, 
tout en permettant un fonctionnement permanent sans pratiquement aucune usure. 


Par rapport aux paliers lubrifiés par liquide, la viscosité bien plus faible des gaz 
réduit encore la perte de puissance, mais exige en même temps des vitesses de 
rotation nettement plus élevées jusqu’à la formation d’un film lubrifiant porteur 
et rend le système plus sensible aux instabilités. Afin de remédier à l'insuffisance 
d'amortissement fluide, une structure constituée de feuilles flexibles introduit du 
frottement sec dans le système et affaiblit ainsi les vibrations auto-entretenues. 
Globalement, cette technologie offre des performances et une fiabilité exceptionnelles, 


1 L’acronyme GFB utilisé ci-après provient du terme anglais Gas Foil Bearing. 
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(a) Dessin technique montrant une VCM d’AiResearch des 
années 1980 avec des GFBs lubrifiés par réfrigérant, utilisées 
dans Lockheed P-3C de l'US Navy (reproduit de [DWHoo)]). 


(b) Photographie d’un GFB open source 
(GFB de type « bump » de 1™ génération) 
conçu par la NASA (adaptée de [NASo7]). 


Figure 1 : Exemples illustratifs de GFBs dans des applications aéronautiques. 


même dans des conditions d'exploitation très exigeantes telles que, par exemple, 
des températures extrêmes. Néanmoins, l'avantage le plus exceptionnel des GFBs 
pour un large éventail d'applications est la possibilité de pouvoir utiliser toute sorte 
de fluide de travail comme lubrifiant, ce qui ouvre la voie à des machines tournantes 
sans huile qui éliminent tout risque d’éventuelle contamination de l’environnement. 


Le développement des GFBs est étroitement lié à leur utilisation fréquente dans 
les systèmes de climatisation des avions, parmi lesquels il faut distinguer entre les 
machines frigorifiques à air (ACMs?) et, comme le montre la Figure 1a, les machines 
frigorifiques à vapeur (VCMs3). Au lieu de l'air ambiant employé dans les ACMs, 
les VCMs utilisent des réfrigérants spéciaux comme l'hydrocarbure halogéné R-245fa, 
permettant de passer par des transitions de phase dans un cycle thermodynamique. 
Comme réalisation possible d’un palier lubrifiable au R-245fa diphasique, la photo- 
graphie de la Figure 1b correspond au design le plus courant des GFBs de type 
« bump » de la 1" génération, que l’on trouve aussi en variantes avec plusieurs lobes. 


L'étude de rotors montés sur des GFBs exige des compétences dans divers domaines, 
notamment la dynamique des rotors, la dynamique des fluides, la thermodynamique, 
la dynamique des structures et la tribologie. Par conséquent, la recherche menée 
à l'échelle mondiale représente une tâche interdisciplinaire qui s'appuie de plus 
en plus sur des outils assistés par ordinateur. Pourtant, l'énorme complexité du 


? L'acronyme ACM utilisé ci-après provient du terme anglais Air Cycle Machine. 
3 L'acronyme VCM utilisé ci-après provient du terme anglais Vapor Cycle Machine. 
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problème rend encore très difficile la mise en place d’une approche de modélisation 
et de résolution à la fois exhaustive et performante. En particulier, la littérature 
manque de modèles minimaux pour les systèmes GFB-rotor lubrifiés par réfrigérant 
qui soient accessibles aux approches monolithiques. Compte tenu de cette lacune 
dans la recherche, le présent mémoire de thèse vise à surmonter les limites de 
faisabilité actuelles afin que, pour la première fois, des études dynamiques de rotors 
puissent contribuer à une compréhension approfondie des VCMs équipées de GFBs. 


Modélisation des systèmes GFB-rotor 


En se basant sur la définition de la géométrie d’un palier et la détermination des 
relations cinématiques fondamentales, le Chapitre 2 établit les préalables nécessaires 
dans l'objectif de pouvoir formuler un problème d'interaction fluide-structure-rotor. 


Le dessin de la Figure 2a montre la section transversale à mi-longueur axiale d’un 
GFB avec trois lobes, dans lequel est monté un des pivots d'arbre du rotor associé. 
Afin d'illustrer les conditions cinématiques pertinentes, la représentation encore plus 
abstraite de la Figure 2b amplifie manifestement les proportions du jeu radial C 
par rapport au rayon R. Le nombre de lobes Np ainsi que leur orientation x, sont 
introduits comme d’autres paramètres géométriques essentiels, créant un modèle qui 
couvre toutes les configurations imaginables de façon aussi universelle que possible. 


La description des champs scalaires ou vectoriels dans l’interstice du palier requiert 
les coordonnées et z, sachant que la dépendance du temps f est aussi considérée. 
Représenté de la manière la plus générale possible, le champ de déformation des 
feuilles est donné par q(ọ, t), tandis que le déplacement d’un pivot peut être exprimé 
par (t) et y(t) ou encore par l’excentricité e(t) combinée avec l'angle de calage T(t). 
La vitesse de rotation du rotor nọ et donc la pulsation wo = 2rrng sont considérées 
constantes, rendant les processus montées—descentes en vitesse quasi-stationnaires. 


Selon une relation trigonométrique appliquée au triangle mis en évidence dans la 
Figure 2b et après linéarisation par rapport à l’excentricité, l'épaisseur du film s'écrit 


Amin = C—e(t) cos[p —I'(t)] — 4(9,#) (1a) 
— [&(t) sin + y(t) cos p] — q(t) (1b) 

a ‚= q (9t) (1c) 

comme superposition du jeu radial, de la contribution sinusoïdale du déplacement 


d’un pivot et du champ de déformation des feuilles. Afin de compléter le couplage 
fluide-structure-rotor, la dérivée temporelle de l’Équation (1) est également formée. 
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Résumé étendu 


(a) Illustration inspirée par le design typique des (b) Géométrie abstraite avec jeu radial largement 
GFBs de type « bump » usuels de la 1 génération. amplifié pour visualiser les relations cinématiques. 


Figure 2 : Dessins en coupe transversale d’un GFB à trois lobes avec un arbre de rotor monté à l’intérieur, 
montrant le fluide (en bleu) qui sépare le pivot tournant (en rouge) de la structure extérieure. Les lobes se 
composent de feuilles supérieures lisses (en vert foncé) supportées par des feuilles ondulées (en vert clair), 
qui se déforment sous l'influence de la pression de l'écoulement et déterminent le contour de l’interstice. 


Enfin, il est possible d'établir un modèle de rotor de type Jeffcott-Laval modifié, selon 
lequel la masse n’est pas intégralement concentrée dans le disque, mais partiellement 
reportée dans les pivots. En tenant compte du poids propre et d’un balourd statique, 
il en ressort un système d'équations différentielles ordinaires qui fait intervenir les 
forces développées par les GFBs en intégrant sur un champ de pression p(@,z,t). 
Grâce à l'approche choisie, on peut également envisager des modèles de rotor plus 
complexes, comme ceux qui deviennent nécessaires pour des conceptions pratiques. 


Théorie de lubrification pour 
les fluides changeant de phase 


La théorie de lubrification classique permet de calculer les distributions de pression 
dans les espaces étroits en simplifiant les équations de la mécanique des fluides. 
Étant donné que l'analyse des GFBs lubrifiés par réfrigérant remet en cause certaines 
des hypothèses sous-jacentes, le Chapitre 3 revoit soigneusement toute cette théorie. 


Si l’on applique les lois de la thermodynamique d'équilibre à un milieu continu, 
tout état possible du lubrifiant est décrit localement de manière unique en spécifiant 
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Théorie de lubrification pour les fluides changeant de phase 


la masse volumique p et la température 0. Ainsi, pour un réfrigérant quelconque, 
toutes les autres variables d’état thermodynamique peuvent étre exprimées sous 
forme d’équations d’état, y compris la pression p, la capacité thermique massique cy, 
le coefficient de compressibilité a, ou encore la viscosité dynamique u. Dans la plage 
de fonctionnement visée, de fortes non-linéarités sont observées pour le R-245fa 
lorsqu’une limite de phase est franchie, ce qui rend inadaptee la loi des gaz parfaits. 


Les profils de vitesse de l'écoulement dans un film mince peuvent être approximés 
avec des équations de Navier-Stokes simplifiées, de manière à ce que leur inclusion 
dans l'équation de continuité donne l'équation de Reynolds pour les gaz non idéaux 


0 wo ð 1 9 [ oh? op d { oh? op 
h) 4 h) = ; 
ot G ) 2 op (e ) R2 a9 e ap dz | 12u dz (2) 
Cette équation aux dérivées partielles décrit directement les variations de pression 


dans le fluide, à condition qu’elle soit tout d’abord fermée par des équations d'état 
thermodynamique et puis couplée à la structure et au rotor en insérant l’Équation (1). 


De plus, étant donné que les propriétés des réfrigérants sont plutöt sensibles a la 
température, il paraît indispensable de simplifier aussi l'équation de l'énergie pour 
un film mince et d’établir une équation de la température pour les gaz non idéaux 


| 4 #098 _ I? (1 ap ae, pa 
Pew) OE” 2 dp 12u \ R20p0@ ` dz dz 


ar 1 (KB ap) of n° op 
PIET? dp =R20¢ \ 12u 9 dz | 12u dz 
u h | 1 /ðp 2 op s 
a ees L | 
Bo( 8 — 9) + Rai Er ah) | (2) . (3) 


Comme dans l'équation de Reynolds, des équations d'état et l'épaisseur du film sont 
incluses, avec le coefficient de transfert thermique fo et la température ambiante Ooo. 


En tant qu’équations aux dérivées partielles fortement non linéaires, l’Équation (2) 
et I’Equation (3) sont d’un caractère mixte parabolique—hyperbolique, ce qui impose 
des restrictions au choix de conditions aux limites pour garder un sens physique. 
Apres avoir identifié tous les domaines dans lesquels des solutions sont a trouver 
ainsi que les bords associés, il faut admettre certaines hypothèses sur l’état du fluide 
autour des GFBs avec pour résultat des conditions aux limites de Dirichlet. Au total, 
il en decoule un modele du fluide sous la forme de plusieurs problemes aux limites, 
dont le nombre précis dépend d’une répartition du domaine en fonction des lobes. 
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Résumé étendu 


Ray > 
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Figure 3 : Dessin du modèle à paramètres concentrés pour un lobe constitué de cinq raidisseurs, réduisant 
la feuille ondulée (« bump foil » en vert clair) à un simple arrangement de tiges, de ressorts et de masses, 
tandis que la feuille supérieure (« top foil » en vert foncé) est interpolée linéairement entre les raidisseurs. 


Modèles performants des feuilles 
avec propriétés de frottement 


Comme la structure à feuilles joue un rôle décisif pour le comportement du système, 
il existe déjà dans la littérature un énorme choix de modélisations correspondantes, 
qui sont dans leur majorité beaucoup trop détaillées et donc assez peu performantes. 


Le dessin de la Figure 3 montre le nouveau modèle à paramètres concentrés présenté 
au Chapitre 4, qui imite les différentes sections de la feuille ondulée en remplaçant 
chacun des raidisseurs par un arrangement de deux tiges rigides et d’un ressort. 
Tout comme les rigidités kg, les masses ponctuelles mg attribuées aux raidisseurs 
sont déterminées sur la base des propriétés de feuilles utilisées dans l’industrie. 
Méme si les fréquences propres des feuilles sont sans intérét quant 4 la dynamique 
du systeme, le maintien de I’inertie s’avere bien utile pour les méthodes numériques. 


Pour la description cinématique du modèle, les coordonnées généralisées x,,1 sont 
introduites comme déplacements des masses ponctuelles, chose particulièrement 
favorable par rapport au frottement à y considérer. De plus, les déplacements yy m 
appartenant aux sommets sont également reliables à cette même représentation 
en exploitant une matrice jacobienne au niveau des vitesses. Entre les raidisseurs, 
la feuille supérieure est simplement approximée par une interpolation linéaire, 
fournissant ainsi un champ de déformation qui délimite l’interstice de lubrification. 


Grâce à la géométrie répétitive et aux nombreuses forces de réaction, les équations du 
mouvement peuvent être facilement dérivées à partir des équations de Lagrange de 
deuxième espèce. Pour chaque lobe, on obtient un système d'équations différentielles 
ordinaires de même forme que pour une simple chaîne d’oscillateurs harmoniques. 
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Analyse numérique des systemes GFB-rotor 


Les forces verticales Fy, correspondent à la sollicitation par le fluide et elles résultent 
d’integrales sur un champ de pression qui sont évaluées dans un voisinage proche. 


Dans les contacts entre les raidisseurs et le fourreau, les forces de frottement avec le 
coefficient ug en adhérence comme en glissement donnent l'expression multivaluée 


{—y8Fnm } si tnm < 0, 
fam € [-uBFnm, +uBFn,m] Si Xn,m = 0, (4) 
{+upF nm} si n,m > 0. 


Sachant que le frottement sous la forme de l’Équation (4) est très difficile à maîtriser 
numériquement, il convient de négliger l’adherence et de la régulariser en utilisant 


x 
fam = UBEnm tanh (=) (5) 

UB 
pour exprimer la loi de Coulomb, ou d’employer un modele de brosse (ou de lame) 
adapté qui soit capable de vraiment reproduire les transitions adherence-glissement. 


Analyse numérique des systemes GFB-rotor 


Avant d’étudier numériquement l’ensemble du systeme, il paraît utile de procéder à 
un adimensionnement cohérent de toutes les équations obtenues de la modélisation. 
Pendant cette procédure apparaissent toute une série de nombres adimensionnels 
qui possedent une signification physique et caractérisent qualitativement le probleme. 


Pour trouver des solutions approximatives de l'équation de Reynolds et de l'équation 
de la température, les dérivées spatiales sont discrétisées au moyen d’une méthode 
des différences finies, un schéma décentré amont étant le plus approprié. Enfin, cela 
transforme le modèle du fluide en un système d’équations différentielles ordinaires. 


La modélisation qui est proposée dans ce mémoire permet de décrire sans ambiguïté 

et à tout moment l’état du système par un seul vecteur d'état fluide-structure-rotor 
T ser et | N 

s=] st st st |encRN, (6) 

appartenant à un espace des phases (2 à grande dimension. Après une simple trans- 

formation des équations précédemment dérivées, on obtient le système dynamique 


D = R{S, Ao}, R:QXxR— R", (7) 
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Résumé étendu 


(a) Solution ponctuelle. (b) Solution périodique. (c) Solution quasi-périodique. 


Figure 4 : Diagrammes de trajectoires transitoires du disque de rotor obtenues par intégration numérique, 
donnant une impression des types de solutions vers lesquels le système dynamique évolue à long terme. 


où la vitesse de rotation adimensionnelle Ag représente le paramètre de bifurcation. 
Ainsi, bien que les sous-modèles pour le fluide, la structure et le rotor soient conçus 
séparément, cette démarche essentielle conduit à une stratégie de solution strictement 
monolithique. Sur la base du Chapitre 5, le code de recherche GFB++ est développé 
en C/C++ pour bénéficier d’une grande vitesse d'exécution et pour profiter du 
paradigme de programmation orientée objet en réalisant une conception modulaire. 


Effets non linéaires dans les systèmes GFB-rotor 


En fonction des paramètres, le système dynamique non linéaire peut se comporter 
très différemment à long terme, ce qui entraîne une classification en différents types 
de solutions selon les graphiques de la Figure 4. La situation est d’autant plus 
compliquée que les solutions peuvent devenir instables et coexistent souvent avec 
d’autres solutions, impossible à détecter par l'intégration temporelle par force brute. 


Une méthode reconnue pour l’investigation systématique des systèmes dynamiques 
est la continuation des solutions, qui est liée aux analyses de stabilité et de bifurcation. 
Pour donner une idée des nombreuses études effectuées, les graphiques de la Figure 5 
montrent l'influence du paramètre d'amortissement adimensionnel DR du rotor sur 
le comportement des systèmes GFB-rotor lubrifiés par réfrigérant en fonction de Ag. 
Évidemment, les scénarios de bifurcation d’après la Figure 5a diffèrent de telle sorte 
que la bifurcation de Hopf (HB) survenue est soit sous-critique, soit supercritique, 
et que celle-ci peut être suivie de jusqu’à deux points selles de cycles limites (LPCs*). 


4 L'acronyme LPC utilisé ci-après provient du terme anglais Limit Point of Cycles. 
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Vitesse de rotation adimensionnelle Ag 


(a) Diagrammes de bifurcation pour les déplacements du disque (en haut) et les fréquences (en bas). 


(b) Trajectoires exemplaires du disque (tracés supérieurs) et des pivots d’arbre (tracés inférieurs). 


Figure 5 : Résultats de la continuation et de l'analyse de bifurcation pour différents amortissements 
(en solide : solutions stables, en pointillé : solutions instables, marqueurs : bifurcations/solutions critiques). 


Résumé étendu 
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(a) Champ de frac. de vapeur (©, = 144, Np = 3). (b) Champ de frac. de vapeur (Op = 144, Np = 1). 


Figure 6 : Diagrammes illustrant la condensation dans des GFBs a un ou a trois lobes. 


Transféré aux applications pratiques, le franchissement d’une HB entraine la perte 
de stabilité d’un point de fonctionnement stationnaire, produisant des vibrations 
auto-entretenues non souhaitables. Tandis que les amplitudes peuvent étre tirées du 
graphique supérieur de la Figure 5a, le graphique inférieur fournit les fréquences 
associées par rapport 4 Ag, la valeur 1% marquant la demi-fréquence de fouettement. 
Enfin, la Figure 5b montre les déplacements statiques et des orbites de cycles limites 
en tant que projections sur les coordonnées du disque de rotor et des pivots d’arbre. 


Une caractéristique particuliére de tous les GFBs lubrifiés par réfrigérant est que des 
transitions de phase entre gaz et liquide peuvent survenir pendant le fonctionnement, 
ce que l’on peut observer pour tous les types de solutions en fonction des paramètres. 
En fonctionnement stationnaire des GFBs à un ou à trois lobes, les graphiques de la 
Figure 6 illustrent la fraction de vapeur dans l’interstice de lubrification sous la forme 
de champs scalaires, le paramètre adimensionnel de gestion thermique ©, = 144 
étant sélectionné de maniere a ce que le fluide ne se réchauffe pas excessivement. 
Comme conséquence de la condensation locale, la distribution de pression change 
considérablement de forme, ce qui finit par nuire a la capacité de charge des GFBs. 


En plus des analyses sur le rôle du frottement sec et sur l'influence d’un balourd, 
quelques études de paramètres sont également réalisées dans le cadre du Chapitre 6, 
dont les principales conclusions sont présentées dans la partie de synthèse ci-après. 


Conclusion et perspective 


Ce mémoire de thése porte sur la dynamique des rotors sur des GFBs lubrifiés par 
réfrigérant et élabore une nouvelle approche de modélisation destinée à permettre 
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Conclusion et perspective 


des analyses de stabilité et de bifurcation. En identifiant les mécanismes d’interaction 
fluide-structure-rotor, il est possible de transformer en un seul systeme dynamique 
trois sous-modeles de complexité plutöt gerable malgré la considération de toutes les 
non-linéarités essentielles. Comme caractéristique principale de la modélisation, 
les équations d'état thermodynamique qui sont directement intégrées dans l'équation 
de Reynolds comme dans l'équation de la température décrivent le comportement 
non idéal d’un réfrigérant typique avec des transitions de phase entre gaz et liquide. 
De plus, le modèle à paramètres concentrés de la structure à feuilles permet la prise 
en compte du frottement sec de plusieurs manières différentes, de la régularisation 
performante de la loi de Coulomb aux modèles de brosse capables de détecter 
les transitions adhérence-glissement. Au total, cela rend accessible le problème 
entier à la théorie mathématique rigoureuse des systèmes dynamiques et permet 
l'implémentation d’un code de recherche monolithique à modules interchangeables. 


Lorsque l’on considère le comportement à long terme des systèmes GFB-rotor 
lubrifiés par réfrigérant, il en résulte l’évolution d'états transitoires vers des solutions 
ponctuelles, périodiques, quasi-périodiques ou chaotiques en fonction des paramètres 
et de la vitesse de rotation systématiquement variée. Pour tous les rotors sans balourd, 
on observe que les positions d'équilibre perdent leur stabilité à une vitesse critique 
à cause d’une bifurcation de Hopf sous-critique ou supercritique, ce qui entraîne 
typiquement l'apparition immédiate de vibrations auto-entretenues non souhaitables. 
Sous les sollicitations dynamiques alors présentes, les GFBs lubrifiés par réfrigérant 
sont particulièrement sensibles aux possibles transitions de phase, un phénomène 
qui peut aboutir à l'effondrement du film lubrifiant et à la défaillance du système. 
En fonctionnement stationnaire, la condensation pose problème aux faibles vitesses, 
bien que les GFBs chargés se réchauffent généralement à tel point que la pression 
de vapeur saturante n’est atteinte nulle part. Afin de rendre les GFBs globalement 
moins sensibles aux oscillations, souvent au détriment de leur capacité de charge, 
le nombre et l'orientation des lobes constituent les propriétés les plus influentes. 
Ensuite, on peut également montrer que la flexibilité et la mobilité des feuilles 
comme des conditions préalables à l'apparition de frottement jouent un rôle décisif 
dans la réduction des amplitudes de vibration jusqu’à atteindre un niveau tolérable. 
Tout à la fin de l'étude, la prise en compte d’un balourd réaliste ne remet pas en 
cause les conclusions précédentes, mais les étend plutôt à des scénarios qui sont 
encore plus complexes car directement liés à des mécanismes d’excitation combinée. 


Enfin, les conclusions du présent mémoire de thèse soulèvent un certain nombre de 
questions de suivi et fournissent ainsi des points de départ pour des travaux futurs, 
portant par exemple sur une amélioration du modèle ou sur les essais expérimentaux. 


Le mémoire de thèse est complété par une liste de symboles et par une bibliographie. 
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1 Introduction 


1.1 Motivation 


Until the present day, the occurrence of detrimental or even destructive vibrations 
remains a serious issue in engineering disciplines that deal with rotating machinery. 
Besides mostly unavoidable synchronous oscillations due to residual rotor unbalance, 
fluid film bearings may additionally induce subsynchronous self-excited vibrations. 
Despite this major drawback, which is known since the early hours of oil lubrication, 
an enormous scientific interest lately develops especially for gas foil bearings (GFBs). 
Based on the same principle of fluid-dynamic lubrication as primitive plain bearings, 
such GFBs are self-acting and thereby rely on the formation of a lubricating wedge 
that builds up pressure due to viscous fluid being dragged into a convergent gap. 
When a fast-spinning rotor shaft is entirely separated from the surrounding structure 
and embedded in a continuous fluid cushion, the absence of solid-to-solid contact 
offers outstandingly low bearing torques and allows for almost wear-free operation. 


Compared to liquid-lubricated bearings, the many times smaller viscosity of gases 
further reduces power losses, but requires significantly higher rotational speeds to 
create a load-carrying fluid film and also makes the system more prone to instabilities. 
As a remedy, the weak fluid damping is supplemented by a compliant foil structure 
that introduces dry friction into the system and thus mitigates self-excited vibrations. 
Altogether, this creates a technology with remarkable performance and reliability 
even under severe operating conditions such as extremely elevated temperatures. 
For the majority of applications, however, the most distinctive advantage of GFBs 
emerges from the usability of practically all kinds of working fluids as a lubricant, 
which ultimately enables oil-free rotating machinery with zero risk of contamination. 


The analysis of rotors on GFBs supposes knowledge covering various fields such as 
rotor dynamics, fluid dynamics, thermodynamics, structural dynamics, or tribology. 
Accordingly, ongoing worldwide research on this topic is an interdisciplinary task 
that currently becomes more and more dependent on predictive computational tools. 
However, the sheer complexity of the problem still poses a considerable obstacle on 
the way toward an all-encompassing but efficient modeling and solution approach. 


1 Introduction 


1.2 State of the Art 


1.2.1 Early Beginnings of GFBs 


Historically, the idea of compliant foil bearings goes back to Blok and Rossum [BR53] 
in the 1950s, even though their concept originally relies on oil as the lubricating fluid 
that builds up pressure between the rotating shaft and a wrapping cellophane foil. 
Some years later, Patel and Cameron [PC57] suggest using a thin steel foil instead, 
which allows for first temperature measurements of such oil-lubricated bearings. 
In the following time, the research works of Baumeister [Bau63] and others at IBM 
establish the connection between foil bearings and gas-lubricated magnetic tapes, 
which at high speeds are separated from read-and-write heads by thin films of air. 
For completeness, it must be mentioned that gas lubrication in itself is not very new 
since it is already addressed by French physicist Gustave-Adolphe Hirn (1815-1890). 
Finally, it is former IBM engineer Gross [Gro63] who has significant influence on the 
pioneering role of US company AMPEX in first developing shaft-supporting GFBs. 
For additional information on the beginnings of compliant foil bearing technology, 
the reader is referred to DellaCorte [Delo8] for a more complete historical outline. 


1.2.2 GFBs in Refrigeration Systems 


Starting from the late 1960s, the commercialization of GFBs is closely related to the 
need for reliable environmental control systems (ECSs) in military and civil aircraft. 
Having to constantly supply crew and passengers with fresh air while maintaining 
the cabin pressure, most jetliners use bleed air that comes directly from the turbines 
and must therefore be cooled down to comfortable temperatures (see, e.g., [BH12]). 
Based on very similar thermodynamic cycles’, there are mainly two variants of such 
cooling packs, namely, air cycle machines (ACMs) and vapor cycle machines (VCMs). 
While VCMs use refrigerants such as Freon or R-245fa that undergo phase changes 
during operation (similar to working principles of domestic A/C and refrigeration), 
common ACMs require nothing but the air itself as a purely gaseous working fluid. 


First practical applications of GFBs in ECSs are reported by Barnett and Silver [BS70] 
from US company Garrett AiResearch, who describe in 1970 the ACM construction 
for the new McDonnell Douglas DC-10 after some tests with retrofitted Boeing 727. 
Moreover, the same authors give important hints on the successful application of 
refrigerant-lubricated GFBs in cryogenic turbocompressors working with Freon. 


t ACMs realize a reverse Joule-Brayton cycle, whereas VCMs realize a reverse Clausius— Rankine cycle. 
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(a) Engineering drawing of AiResearch VCM from the 1980s (b) Photography of an open-source GFB 
on refrigerant-lubricated GFBs, which are first installed on (bump-type GFB of the first generation) 
Lockheed P-3C of US Navy (reproduced from [D WH9ol). made by NASA (adapted from [NASo7]). 


Figure 1.1: Illustrative examples for GFBs in aeronautical applications. 


During the 1970s, the leaf-type GFBs invented and patented by Garrett AiResearch 
are more and more competed by bump-type Hydresil GFBs first developed by MTI 
and later employed by Hamilton Standard in ACMs of aircraft such as Boeing 747. 
In this connection, the report from 1978 by Ruscitto, McCormick, and Gray [RMG78] 
represents certainly the most comprehensive study that is available from this period. 
Over the years, both leaf-type and bump-type GFBs become firmly established in 
the aviation industry as can be seen from a retrospective article by Agrawal [Agrg7]. 
In a wider context, GFBs are slowly gaining ground in other oil-free applications 
that notably include turbochargers or combustion turbines [Howgg, HBRıo, Del12]. 


Regardless of the bearing type, manufacturers of large aircraft traditionally prefer 
ACMs over VCMs because of their particularly simple construction (see, e.g., [MSo8)). 
However, critical reviews of this situation by Dexter, Watts, and Haskin [DWHogo] 
in 1990 and by Linnett and Crabtree [LC93] in 1993 outline advantages and future 
opportunities of VCMs especially in combination with refrigerant-lubricated GFBs. 
According to them, ongoing trends toward bleedless turbines and electric aircraft 
can make VCMs the more attractive, more efficient, and more cost-effective solution 
since the primary reason for using ACMs is eliminated together with the bleed air. 
In contrast to less promising attempts with pneumatically driven VCMs in the past, 
the availability of powerful electric motors nowadays allows for the construction of 
such systems that operate reliably even during hot ground phases (see, e.g., [BH12]). 
One step in this direction is documented by the engineering drawing in Fig. 1.1a, 
which shows an electrically driven VCM running on refrigerant-lubricated GFBs 
constructed by AiResearch for military applications and installed on Lockheed P-3C. 


1 Introduction 


Obviously, the reduction of masses and volumes plays a decisive role in aeronautics, 
which typically leaves little flexibility to engineers when designing ACMs or VCMs. 
With GFBs, this usually means that the working fluid is also used for lubrication 
since additional lubricant circuits with the required sealing may then be economized. 
Especially in VCMs, refrigerant lubrication even seems without alternative because 
the slightest contamination with oil is suspect to drastically reduce thermodynamic 
performance according to general investigations by Youbi-Idrissi and Bonjour [YBo8]. 
As of today, large commercial aircraft are still equipped with ACMs as standard, 
including the extensively electrified and bleedlessly designed Boeing 787 Dreamliner. 
Nevertheless, the ongoing research in leading companies such as Liebherr-Aerospace 
is documented in a report by Faleiro, Herzog, Schievelbusch, and Seung [Fal+o4] 
and in the more GFB-specific work of Garcia, Bou-Said, Rocchi, and Grau [Gar+13]. 
Elsewhere, efforts in this direction are also vaguely suggested but not made available. 


The selection of refrigerants for industrial applications is subject to many constraints, 
which are outlined in the detailed overview article by Calm and Didion [CD98]. 
Besides the fundamental issues of toxicity and flammability, environmental concerns 
including the depletion of the stratospheric ozone layer seem directly related to the 
use of chlorofluorocarbons or hydrochlorofluorocarbons such as Freon refrigerants. 
Following the Montreal Protocol, most countries substitute chlorinated refrigerants 
by alternatives such as hydrofluorocarbons, which cause no harm to ozone molecules. 
For the specific application in aeronautical VCMs, the hydrofluorocarbon R-245fa 
(1,1,1,3,3-pentafluoropropane C,H;F-) is found a particularly promising candidate 
from the thermodynamic point of view according to the thesis of Röyttä [Röyog]. 
Moreover, other studies on refrigerant-lubricated GFBs in VCMs by Garcia [Gar12] 
and Bouchehit [Bou17] are principally conducted with substance data of R-245fa. 
Recently, however, hydrofluorocarbons also fall into disrepute due to them being 
greenhouse gases that are suspect to contribute to human-induced climate change. 
Future replacements are often seen in hydrofluoroolefins or hydrochlorofluoroolefins 
as investigated in a study by Eyerer, Dawo, Kaindl, Wieland, and Spliethoff [Eye+19]. 


The photography in Fig. 1.1b shows a typical bump-type GFB of the first generation 
that corresponds to the NASA’s open-source design by DellaCorte, Radil, Bruckner, 
and Howard [Del+08], which is created after the expiry of patents to make GFBs 
accessible to public research. Besides this classical construction, multi-pad layouts 
are focused on from the very beginning by Heshmat, Shapiro, and Gray [HSG82] 
and until today in studies by Osmanski, Larsen, and Santos [OLS17] and elsewhere. 
Other recent bump-type GFBs are mostly credited to the specialized company MiTi, 
where Heshmat, Walton, and Cördova [HWC18] lately mention a sixth generation. 
Finally, there are some other patented solutions with alternative designs such as the 
wing-type GFBs by Swanson and O’Meara [5018], which are not further considered. 


1.2 State of the Art 


1.2.3 Analysis of Rotors on GFBs 


Since the earliest studies on GFBs, numerical analyses and simulations prove to be 
indispensable tools for the design and evaluation of such systems, not only because 
experiments are expensive, complex, and inflexible, but also because many of the 
interesting quantities in fluid and foil structure seem hardly or not at all measurable. 
The following gives a brief overview of developments in this field that are relevant 
for the topic of the present thesis, whereas many additional aspects are covered by 
the systematic review recently written by Samanta, Murmu, and Khonsari [SMK19]. 


An important foundation is laid with a book by Walowit and Anno [WA75] in 1975, 
which precedes a much acclaimed paper of Heshmat, Walowit, and Pinkus [HWP83]. 
Their description of the fluid goes back to the work of Reynolds [Rey86] from 1886, 
whose famous equation is extended by Harrison [Har13] to include compressibility. 
Besides the discussion of suitable boundary conditions for the Reynolds equation, 
the mentioned publications on GFBs propose a representation of the foil structure 
using a simple elastic foundation? model calibrated with analytical bump stiffnesses. 
Due to their simplicity and efficiency, foundation models are popular until today, 
in particular after undergoing some improvements and refinements over the years. 
Most importantly, the addition of viscous damping by Peng and Carpino [PC93] 
allows for a basic consideration of frictional energy dissipation mechanisms with the 
equivalent coefficients from Peng and Carpino [PC94] or Ku and Heshmat [KH94]. 
Another class of amendments concerns definitions of the foundation stiffness itself, 
which could vary along the circumferential direction in the spirit of lordanoff [lor99], 
which could be averaged in the axial direction as done by Peng and Khonsari [PKo4], 
or which could obey to a progressive spring characteristic as suggested by Baum, 
Hetzler, Schröders, Leister, and Seemann [Bau+21] for time-efficient computations. 


During the last 15 years, the increasing performance of available computers correlates 
with the appearance of more advanced foil structure models that include dry friction. 
While there are interesting attempts with three-dimensional finite-element models 
by researchers such as Lee, Kim, and Kim [LKKo8] or San Andres and Kim [SKog], 
the majority of investigations make use of two-dimensional lumped-element models. 
Pioneering work in this context is attributed to Le Lez, Arghir, and Fréne [LAFo7], 
who set up a network of linked springs to represent interactions within the structure. 
Simplifying the structural geometry even further, Feng and Guo [FG14] present the 
inclusion of friction models that overcome limitations of equivalent viscous damping. 
Other promising strategies of reproducing friction characteristics are reported by 
Barzem, Bou-Said, Rocchi, and Grau [Bar+13] or Larsen, Varela, and Santos [LVS14]. 


? Especially in the German language literature, this idealization is also referred to as a Winkler foundation. 


1 Introduction 


Finally, there is another sophisticated model by Arghir and Benchekroun [AB19], 
which takes into account foil detachment based on a contact mechanics approach. 


Mostly parallel to the evolution of foil structure models, major efforts are also made 
by several research groups with regard to the thermal modeling of the fluid behavior. 
While there are already complete thermohydrodynamic descriptions of oil-lubricated 
bearings in the 1980s given by Boncompain, Fillon, and Fréne [BFF86] and others, 
practically all GFB models during this period are based on isothermal assumptions. 
It is only at the turn of the millennium that Salehi, Swanson, and Heshmat [SSHo1] 
publish one of the first comprehensive studies solving an energy equation for GFBs, 
but under assumptions that decouple the temperature from the Reynolds equation. 
Subsequently, mainly two approaches are pursued within fully coupled analyses, 
using either full three-dimensional energy equations like Peng and Khonsari [PKo6] 
or cross-film averaged temperature models as done by San Andrés and Kim [SK10]. 
One of the latest trends is marked by the promising investigations of Garcia [Gar12], 
Kim [Kim15], or Bouchehit [Bou17], who add constitutive equations for non-ideal 
s-CO, or R-245fa, which even allows to reproduce phase transitions in the latter case. 
Another thorough discussion of vapor-liquid phase transitions in R-245fa is given by 
Guenat and Schiffmann [GS19] regarding applications under isothermal conditions. 
Furthermore, there is growing interest in thermal fluid-structure-rotor interaction 
to understand runaways and instabilities described by Samanta and Khonsari [SK18]. 
An experimental investigation with a particular regard to rotor temperatures during 
run-up and operation is reported by Mahner, Bauer, Lehn, and Schweizer [Mah+19], 
while measured foil temperatures are documented by Zywica and Bagiriski [ZB19]. 


Even though designers of rotating machinery are in particular interested in the 
dynamic performance of a rotor, surprisingly few publications focus on this topic. 
Instead, most of the aforementioned studies assume steady-state operating conditions 
and make use of iterative solution strategies to identify converged system states. 
Then, it is often attempted to characterize rotor vibrations using linearized stiffness 
and damping coefficients obtained from the perturbation method of Lund [Lun68]. 
However, Bou-Saïd, Grau, and Iordanoff [BGIo8] can show that linearized approaches 
possess a very limited range of validity and overlook important nonlinear effects. 
Calculating a large amount of rotor trajectories through numerical time integration, 
Bhore and Darpe [BD13] give systematic insights into the rich nonlinear dynamics of 
GFB-rotor systems by consideration of orbit plots, FFT spectra, and Poincaré maps. 
Finally, the feasibility of bringing all governing equations into the general form of a 
mathematical dynamical system is first shown in 2014 by Bonello and Pham [BP 14]. 
Similar approaches are employed in recent articles by Gu, Ma, and Ren [GMR17], 
Leister, Seemann, and Bou-Saïd [LSB19], or Osmanski, Larsen, and Santos [OLS2o0], 
focusing on transient responses, bifurcation behavior, or stability limits, respectively. 


1.3 Objectives and Outline 


Notwithstanding the many advantages of such dynamical system formulations, 
research by Osmanski, Larsen, and Santos [OLS18] or Pronobis and Liebich [PL19] 
also aims to improve perturbation methods for more efficient parameter studies. 
Avoiding classical perturbation entirely, the latest publication of Bonello [Bon2o] 
applies a new method for the efficient generation of numerical Campbell diagrams. 
Often there are also conflicting goals between the desired model complexity and the 
transferability into a dynamical system, which makes Bou-Said, Lahmar, Mouassa, 
and Bouchehit [Bou+20] conclude that monolithic approaches are not always feasible. 


1.3 Objectives and Outline 


In the ubiquitous field of tension between the complexity of multiphysics models 
and the request for applicability of strongly coupled monolithic solution approaches, 
the present thesis is meant to explore and push forward existing limits of feasibility. 
To this end, a first goal consists in the proper identification of interaction mechanisms 
in all directions between the three involved domains of fluid, foil structure, and rotor, 
giving the most general framework for a unified description of the overall system. 
As another goal, the conducted derivation of a modified Jeffcott-Laval rotor model 
not only provides a foundation for subsequent investigations but is also intended 
to generally establish the link to classical rotor dynamics with nonlinear bearings. 
Most importantly, considering the particular case of refrigerant-lubricated GFBs, 
the primary modeling objective consists in finding a suitable compromise to fill the 
currently prevailing gap between overly complex and overly simplistic approaches. 
Taking into account non-ideal gas behavior including vapor-liquid phase transitions, 
a novel thermo-gasdynamic fluid model is to be presented after rigorous derivation 
of a modified Reynolds equation combined with a thin-film temperature equation. 
Also, the proposed lumped-element foil structure model aims to satisfy a demand 
for reproducing dry friction and bump-bump interaction while remaining superior 
in computational efficiency compared to most other attempts featuring such details. 
Altogether, the presented work is original in showing how an extensively modeled 
refrigerant-lubricated GFB-rotor system can be represented as a dynamical system, 
making it possible for the first time to perform stability and bifurcation analyses. 
Following the development of a numerical research code as one part of the project, 
several effects and phenomena serve as the starting points for closer investigations 
and thereby lead to novel findings especially with regard to the nonlinear dynamics. 
While previous studies on refrigerant lubrication just consider exemplary equilibria, 
this completes the state-of-the-art research by systematic performance overviews and 
is supposed to deliver a better understanding of how to prevent or reduce vibrations. 


1 Introduction 


This thesis is articulated in seven individual chapters including the present Chapter 1. 


Chapter 2 provides a generalist introduction to the modeling of GFB-rotor systems. 
Defining the bearing geometry and identifying fundamental kinematic relationships, 
this creates prerequisites for setting up a fluid-structure-rotor interaction problem. 
Referring to the standard modeling approach of most frequent use in rotor dynamics, 
a modified Jeffcott-Laval rotor model that accounts for GFB forces is then considered. 


Chapter 3 revisits the theory of fluid film lubrication for phase-changing refrigerants. 
Having established a thermodynamic framework with suitable equations of state, 
the derivation of a modified transient Reynolds equation is demonstrated in detail. 
Confronted with the temperature sensitivity of phase transitions, an energy equation 
completes the boundary value problems to be solved for given boundary conditions. 


Chapter 4 presents a lumped-element foil structure model that considers dry friction. 
Essentially, the reduction of the continuum to a simple spring—mass arrangement 
allows for an elegant formulation of kinematic relationships and equations of motion. 
Moreover, it is shown how dry friction can be included in varying levels of complexity, 
reaching from regularized Coulomb friction models to elasto-plastic bristle models. 


Chapter 5 is dedicated to the efficient computational analysis of GFB-rotor systems. 
Besides the definition of a nondimensional problem formulation, the solution strategy 
encompasses a discretization of fluid equations by upwind finite-difference schemes. 
Altogether, the fluid-structure-rotor submodels form a single dynamical system, 
for the treatment of which a tailor-made modular research code GFB++ is developed. 


Chapter 6 guides through the most important nonlinear effects and phenomena. 
After an illustrative introduction into the mathematical theory of dynamical systems, 
the description and discussion of relevant findings is divided into three sections. 
Firstly, the particular features of refrigerant lubrication with phase transitions are 
focused on with regard to both steady-state operation and vibrational self-excitation. 
Secondly, the influence of some foil structure design parameters and the importance 
of dry friction are thoroughly investigated using stability and bifurcation analyses. 
Finally, the addition of rotor unbalance gives an outlook to quasi-periodic behavior. 


Chapter 7 concludes the thesis with a summary and a perspective to future research. 


2 Modeling of GFB-Rotor Systems 


In the field of theoretical rotor dynamics, it is well-known that assumptions regarding 
the bearings have a crucial influence on the predicted behavior of the overall system. 
To create a basic understanding of the modeling of GFB-rotor systems in particular, 
this chapter begins with a general overview of the bearing geometry and establishes 
fundamental kinematic relationships between rotor, fluid film, and foil structure. 
Understanding the relevant interaction mechanisms, it is finally possible to consider 
a modified Jeffcott-Laval rotor model that accounts for state-dependent GFB forces. 


2.1 Geometry of Gas Foil Bearings 


The schematic drawing in Fig. 2.1a shows a cross-section through the circumferential 
centerline of a typical bump-type GFB together with the mounted rotor shaft journal. 
In the considered case of a simple multi-pad bearing similar to those fabricated 
by Heshmat et al. [HSG82], the compliant metal foil structure is split into several 
identical segments, each of which consists of a smooth top foil layer resting on a 
uniformly corrugated bump foil strip. To allow for some relative motion within 
the structure while keeping the pads in position, these foils are loosely inserted 
into the bearing and spot-welded to the sleeve only at their trailing edges. For the 
sake of generality, design modifications such as axially varying structural profiles or 
irregularly corrugated bump foils are not of interest within the scope of this thesis. 


As a matter of principle, the load-carrying fluid film forming between the journal and 
the top foil during operation is several orders of magnitude thinner than the nominal 
dimensions of the bearing [Agrg7]. For the illustration of kinematic relations within 
this hardly perceivable lubrication gap, the sketch of the cross-sectional geometry 
depicted in Fig. 2.1b greatly exaggerates the proportions of the fluid film and the 
magnitude of structural deformations. In addition, the graphical representation of 
bump foils and top foils is abstracted to an extent that the depicted multi-pad layout 
can easily be transferred to GFBs with other configurations of bumps and pads, 
including conventional single-pad bearings and limiting cases of full-arc bearings. 
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(a) Illustrative drawing inspired by typical designs (b) Abstracted geometry with greatly exaggerated 
of widespread first-generation bump-type GFBs. lubrication gap to see relevant kinematic relations. 


Figure 2.1: Sketches of cross-sections through an exemplary three-pad GFB with mounted rotor shaft, 
showing how the fluid film (in blue) separates the rotating journal (in red) from the surrounding structure. 
The pads consist of smooth top foils (in dark green) supported by corrugated bump foils (in light green), 
which deform under the influence of flow-induced pressure and affect the shape of the lubrication gap. 


The bearing sleeve is modeled as a rigid hollow cylinder with the height to be 
imagined perpendicular to the drawing plane and equal to the axial bearing length L. 
Regardless of the investigated rotor model, the shaft journal with radius R is then 
represented by a rigid cylinder, which is assumed to remain perfectly aligned with the 
bearing axis. According to a numerical study conducted by Carpino et al. [CPMo4], 
slight journal misalignment is compensated by adaptations of the foil structure and 
rubbing at the axial edges of the bearing is thus prevented. Consequently, provided 
an axially coupled structural model with deforming edges, numerical predictions 
are not expected to differ qualitatively whether skew is explicitly considered or not. 


Given that the inner diameter of a GFB prior to mounting is designed to be slightly 
smaller than the nominal diameter of the journal, the foil structure is preloaded 
upon assembly of the rotor depending on the actual interference. This makes the 
existence of an initial clearance less obvious than in the case of rigid bearings 
because the top foil is not clearly separated from the journal when the machine is 
at rest. Nevertheless, a minimum amount of play always remains in the system 
due to the foils not being completely seated against the sleeve and not perfectly 
conforming to the circular shape of the bearing [RHDo2]. In practical applications, 
this clearance can only roughly be estimated from load-deflection curves obtained 
with a non-rotating dummy shaft pushing the foils outward, where a sudden rise 
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in stiffness indicates the available space to be exceeded [RMG78]. To make matters 
worse, such static measurements must be considered with caution since rotor shaft 
growth due to centrifugal forces and thermal expansion may aggravate the issue of a 
hardly quantifiable clearance [KS06]. For modeling purposes, on the other hand, the 
commonly adopted approach consists in the introduction of a well-defined radial 
clearance C > 0 as an empirical representation of the effective assembly preload. 
In the literature, values on a length scale of tens of micrometers are suggested for 
this parameter, whereas shaft journal and bearing dimensions R and L are typically 
found on a length scale lọ of several centimeters [Hes94]. This finding introduces a 
small geometric parameter 


e=— <1 (2.1) 


into the problem, which is expected in the characteristic order of magnitude of 1073 
for the two given length scales and justifies GFBs to be treated within the framework 
of fluid film lubrication theory as established by Reynolds [Rey86] and others. 


As an obvious consequence of the foil structure being arranged in a circular pattern, 
the number of pads Np in a typical bump-type GFB is always equal to the number of 
fixation gaps at which the individual foil segments are welded to the bearing sleeve. 
While the drawing in Fig. 2.1a illustrates this fact for an exemplary three-pad bearing 
with Np = 3 and as many fixation gaps, such an equivalence of pads and gaps holds 
even for designs without gaps if these are considered zero-pad bearings with Np = 0. 
In the simplified model according to the sketch in Fig. 2.1b, the abstracted fixation 
gaps for any Np > 1 have no spatial extension and are evenly distributed around the 
circumference, which means they are periodically spaced by the pad sector angle 


_ 20 


Ax = N (2.2) 


To complete this most universal description of a multi-pad GFB, the orientation of an 
arbitrarily chosen first fixation gap is prescribed by the angle x1, which implies that 


Xn =Xit(n-1)Axy, n=1,...,Np (2.3) 


denotes the respective angular position of the n-th gap in relation to the vertical axis. 
In the considered multi-pad GFBs, each of the Np > 1 bump foil strips consists of 
the same number Ng of identical bumps, which results in a total of NpNg bumps. 
Based on the pad sector angle Ax from Eq. (2.2) as shown in the sketch in Fig. 2.1b, 
the circumferential extension of the undeformed bumps can therefore be described 
by the bump sector angle 


Ay = Ne = NN (2.4) 
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Shell 


Figure 2.2: Sketch of the circumferentially “unwrapped” lubrication gap for an exemplary three-pad GFB, 
which is strongly distorted to illustrate coordinates and contributions to the effective fluid film thickness. 


Moreover, Eq. (2.3) giving the angular position Xn of the leading edge of the n-th pad, 
it is possible with Eq. (2.4) to specify the angular position of the associated m-th bump 


Pam =Xn+(m—1)Ap, n=1,...,Np, m=1,...,Np. (2.5) 


Since Eq. (2.5) obviously refers to the midpoint of the preceding connecting section, 
the apex of the undeformed m-th bump on the n-th pad is located at Py + AY /2. 


2.2 Coordinates and Kinematics 


For the physical description of a GFB in time (denoted by t) and space, the bearing 
sleeve represents not only the system boundary but also an inertial frame of reference. 
In view of the small clearance-to-radius ratio C/R < 1 that is implied by Eq. (2.1), 
the center of curvature of the lubrication gap appears sufficiently distant on the 
length scale of film thickness to henceforth neglect this curvature (see, e.g., [Sze1o]). 
When mentally “unwrapping” the fluid film as illustrated by the sketch in Fig. 2.2, 
the domain can be represented by a set of Cartesian coordinates (x,y,z) associated 
with the origin F and an orthonormal basis {ex, ey, ez}, which are fixed to the inertial 
frame of reference. In doing so, the circumferential coordinate x is identified with an 
arc length that is directly related to the angular coordinate p by the transformation 


x = Roy. (2.6) 


On account of Eq. (2.6), field quantities in the lubrication gap such as pressure can 
be formulated as functions of pọ while being evaluated using Cartesian coordinates. 


Getting back to the sketch in Fig. 2.1b, it follows immediately from the idealized 
assumption of an initial clearance C around the journal of radius R that the inner 
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boundary of the undeformed foil structure is characterized by the radius R + C. 
When considering GFBs under operating conditions, it is known from experimental 
observations by Ruscitto et al. [RMG78] that the entire foil structure deforms almost 
uniformly along the bearing axis even in the case of significant fluid pressure 
gradients in the same direction. Inspired by this practical evidence of high axial 
stiffnesses, a comparative numerical study by San Andrés et al. [SKog] confirms 
that the foils of such bearings can be represented with sufficient accuracy by plane 
structural models under axially averaged loading. Hence, the generic deformation 
field q(,t) is introduced as being a priori independent of the z-coordinate and 
describes, without discussing in detail the kinematics of foil structure models for the 
time being, the effective shape of the top foil as the outer boundary of the fluid film. 


For a kinematic description of the rotor, in which this “unwrapped” view on the 
bearing is not relevant, another orthonormal basis {e;, ey, ez} is fixed to the inertial 
frame of reference as shown by Fig. 2.1b in the center of the sketch. Furthermore, 
the rotation by an angle of m — T(t) around the positive eç-axis transforms this 
system of unit vectors into a time-dependent orthonormal basis {e,(t),er(t), —ez} 
that follows the journal motion. In this way, referring to the origin © located at the 
center of the bearing, the planar position vector 


r(t) = C(t)ee + (ten (2.7a) 
= e(t)ee(t) = e(t) |- sin T (t)eg — cos T (t)ey] (2.7b) 


for the center of the journal M is determined either by the Cartesian coordinates &(t) 
and (+) as in Eq. (2.7a) or by the polar coordinates e(t) and T(t) as Eq. (2.7b) shows. 
While the transformation rules from polar coordinates to Cartesian coordinates 


&(t) = —e(t) sin T (t), y(t) = —e(t)cosT (t) (2.8) 


are obtained by simply equating the coefficients in Eqs. (2.7a) and (2.7b), the inverse 
transformation is slightly more delicate. It consists in a geometrical identification of 
the eccentricity 


e(t) = ee? + nn)? (2.9) 


with the Pythagorean theorem and, defining a two-argument trigonometric function 
arctan2: IR*\ {(0,0)} — (—7, x] to distinguish the quadrants, of the attitude angle 


arccos Pen -m if é(t)>0, 
= ir 5 = 2 if ¢(t) =0 
T(t) = arctan2[—¢(t), -—4(t)] = undefined al (2.10) 
TT — arccos os otherwise. 
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Since the eccentricity and the attitude angle provide a very natural way of describing 
the displacement of the shaft journal relative to the bearing, polar coordinates are 
often the first choice in the discussion and evaluation of results. On the other hand, 
it is already obvious from Eq. (2.10) that polar coordinates bear the disadvantage 
of not being properly defined at their pole O, which makes Cartesian coordinates 
the better choice for the actual rotordynamic calculations [Bau18]. Finally, supposed 
that the constant rotational speed of the rotor no is given in terms of SI units such 
as 1/s rather than in rpm, the angular velocity vector writes 


w = —woez = —271N9ez, (2.11) 


with wo = 2rıng being subsequently referred to as the angular frequency of the rotor. 
When reproducing machine run-ups and coast-downs in spite of wg being constant, 
such transient procedures must be slow enough to be considered quasi-stationary. 


2.3 Fluid-Structure-Rotor Interaction 


When applying the law of cosines to the triangle OMB highlighted in the sketch in 
Fig. 2.1b, where B denotes an arbitrary point at the inner boundary of the fluid film, 
the immediately emerging equation 


R? = e(t)? + [R +C -hlg t) -ql D]? 
— 2e(t)[R+C —h(9,t) —q(9,t)] cos[p-T(t)] (2.12) 


is quadratic in the fluid film thickness field h(p,t) and possesses only one physically 
meaningful solution. Developing the corresponding Taylor series with respect to the 
eccentricity e(t), which is intuitively found on the same length scale as C in Eq. (2.1), 
all terms of orders greater or equal to [e(t)/R]? appear entirely negligible compared 
to unity (see, e.g., [Sze10]). Finally, the linearized fluid film thickness field writes 


h(g,t) =C —e(t) cos[p —T(t)] — (pt) (2.13a) 
=C- [(t) sing + y(t) cos p] — q(p,t) (2.13b) 
= h:(p,t) — q(p;t) (2.130) 


and constitutes a superposition of the nominal clearance, the sinusoidal contribution 
of the displaced rotor shaft journal, and the deformation field of the foil structure. 
To abbreviate the polar representation of the shaft journal position in Eq. (2.13a) or 
the alternative Cartesian representation in Eq. (2.13b), which are directly related by 
trigonometric identities, the most general form in Eq. (2.13c) introduces the film 
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thickness field h, (ọ, t) of a virtually rigid bearing. Looking at the sketch in Fig. 2.2, 
the radial extension of the fluid domain is then given in terms of the y-coordinate by 


h(p,t) Sy < —q(p;t) (2.14) 


since the origin F of the fluid coordinate system is located right on the contour line 
of the undeformed foil structure. In direct analogy to the film thickness field itself, 
the corresponding time derivative field 


ð 


111122) = — let!) cos [ep — T(#)] + e(t)F(t) sin[p — ro] — ag.) (2.15a) 
= [sin + (5) cos 9] — Zalo») (2.156) 
= Žr (g,t) - Sa(g.t) (2150) 


can also be expressed, depending on the context, in any of the discussed forms. 
Unsurprisingly, it is noted that all film thickness relations derived in this paragraph 
are independent of the z-coordinate under the aforementioned assumptions of axially 
uniform foil structure deformations and negligible shaft journal misalignments. 


In contrast to what may be predicted by such an idealized model, the lubricant film 
of a properly operating GFB cannot become arbitrarily thin in practical applications. 
Instead, depending directly on the average surface roughnesses, a certain distance ho 
must always be maintained between the foil structure and the rotor shaft journal in 
order to avoid wear and damage [DVoo]. Consequently, having introduced 

hmin(t) = min h(ọ, t) (2.16) 

vera] 

to characterize the remaining film thickness at the current point of closest approach, 
the relevance of calculated results may be checked at any time t against the condition 


hmin(t) > ho. (2.17) 


Tribologically speaking, the statement of Eqs. (2.16) and (2.17) corresponds to a 
demand for full fluid film lubrication, in which asperity contacts play no role at all. 


The fluid acting as a link between foil structure and rotor, the pressure field p(q, z, t) 
that builds up depends primarily on structural deformations and deformation rates 
as well as on shaft journal displacements and displacement velocities. Consequently, 
the governing equations of the lubricant can be formulated elegantly in terms of 
Eqs. (2.13) and (2.15), which consolidate all those coupling quantities in two fields. 
Conversely, the pressure field is identified as the immediate cause of foil deflections 
on the one hand and of bearing forces acting upon the rotor shaft on the other hand, 
which finally results in a fully coupled fluid-structure-rotor interaction problem. 
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Top foil 
Bump foil 


Figure 2.3: Sketch of the modified Jeffcott-Laval rotor model running horizontally in symmetrical GFBs, 
for which it is advantageous to introduce additional mass particles (in dark red) for the two shaft journals. 


2.4 Modified Jeffcott-Laval Rotor Model 


In order to create a basic understanding of how the peculiarities of GFBs may affect 
the dynamics of any rotor, it is reasonable to consider a rather unspecific rotor model 
that shows well-known behavior in combination with simpler bearing descriptions. 
Here, as in most of the relevant literature on rotor dynamics (see, e.g., [GNPo2)), 
the classical Jeffcott-Laval rotor model represents an appropriate reference for such 
fundamental investigations and it reveals a good compromise in terms of complexity. 
Due to the modular modeling approach, it is almost trivial to complicate or to further 
simplify the rotor model as required for any conceivable application at a later time. 


As can be seen from the schematic sketch in Fig. 2.3, the considered rotor model 
consists of an unbalanced disk that is mounted in the middle of a linearly elastic shaft. 
In this lumped-element approach, the parameter kr characterizes the shaft stiffness, 
while dr represents some non-rotating viscous damping that the surrounding fluid 
opposes lateral disk movements and mp describes the total mass of the entire rotor. 
Here, unlike the original Jeffcott-Laval rotor model, a small portion of mass Ermr/2 
is shifted to either of the shaft journals with the mass distribution parameter ER < 1, 
so that only the remaining mass (1 — ER)"R can actually be attributed to the disk. 
Adding to the physical justification of a model that also accounts for the shaft mass, 
a positive side effect is the resulting system of ordinary differential equations (ODEs) 
being much more convenient than differential-algebraic systems of equations (DAEs). 
The only external load of the horizontal rotor is the weight of the disk and of the 
shaft journals with the gravitational acceleration g acting in the negative e,-direction. 
Finally, giving probably the most radical simplification compared to real applications, 
the rotor is assumed to remain perfectly symmetrical with respect to the two GFBs, 
which ignores conical mode shapes and also gyroscopic effects from the beginning. 
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Figure 2.4: Free-body diagrams of the disk (on the left) and of both symmetrical journals (on the right) 
expressed as one, from which the equations of motion for symmetric mode shapes can readily be derived. 


The sketch in Fig. 2.4 shows free-body diagrams of the shaft journals and of the disk, 
but also clarifies the introduced rotor coordinates and kinematic relations in general. 
Just as Eq. (2.7a) describes the position of a journal center in Cartesian coordinates, 
the position of the center of the disk M’ is specified by &r(t) and yr(t), which gives 
a planar position vector 

rr(t) = Gx(teg + mr (Hey. (2.18) 
In the presence of a static unbalance, for which the center of mass M! is displaced 
from M’ by some eccentricity er, the former is located by another position vector 


rR (t) = Čr (t)ez + Nr,» (t)ey- (2.19) 


If t = 0 corresponds to the instant at which M! is found just vertically below M’, 
the coordinates of the center of mass in Eq. (2.19) are related to those in Eq. (2.18) by 


ČR (t) = p(t) — er sin(wot), NR, (t) = R(t) — er cos(wot) (2.20) 
and depend explicitly on time using the constant angular frequency wo of the rotor. 
With the rotational motion being prescribed according to Eq. (2.11), the rotor model 
has four remaining degrees of freedom that can be identified with ¢, 7, &r, and yp, 
while the bearing force components Fz and F} can be doubled due to symmetry. 


Expressing the total acceleration of the center of mass by deriving Eq. (2.20) twice, 
Newton’s second law of motion applied to the free-body diagrams in Fig. 2.4 yields 


Epmpé + kr (č — CR) = 2Fe, (2.21a) 

ErmRÄ + kr (mM — yr) = 2F, — ERMRg, (2.21b) 
(1 — Ep) még + drép + ke (Ex — €) = —(1 — Sp) menu sin (wot), (2.21¢) 
(1 — Er) mrifp + drig + kr (MR — 7) = —(1 — Er)mrerwg cos(wot) 


— (1- Er)mrg.  (2.21d) 
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When finally reordering Eq. (2.21), it is convenient to introduce the two parameters 


(Re AR 
WR = imp’ Re (2.22) 


usually referred to as natural angular frequency and damping decay rate, respectively. 
For the sake of generality, however, it is not intended to derive numerical values for 
these two parameters from the exact design of any real-world rotating machinery. 
Then, the system of equations of motion for the modified Jeffcott-Laval rotor writes 


gl el 
6 0 ER —1 0 1-5 0 G 
il, 2p lol oR oo | | 9 
ërR| 1-Ep |éx| 1-ER| -1 0 1 0 ER 
ÜR ÜR 0 —1 0 1 AR 
le | at 0 
= 7 2 a 

Ermg | 0 RED sin (wot) 8 lol: (2.23) 

0 cos(wot) 1 


While the left-hand side of Eq. (2.23) describes coupled harmonic oscillators that are 
influenced by the mass distribution parameter Er, the terms on the right-hand side 
are, in that order, nonlinear bearing forces, unbalance forcing, and gravitational load. 
With Eq. (2.23), the limiting case of a rigid rotor with degrees of freedom & and y is 


ë ë| 2 Fe) sin(wot)|  |0 
1 OMR | enw aaa 8 a 2 (2.24) 


~ | + 2dr 
which is a useful simplification for rotor shafts that are much stiffer than the GFBs. 


The bearing force that results from a pressure field p around the shaft journal can be 
obtained by simply integrating elementary pressure forces over that entire surface. 
For a sufficiently thin film according to Eq. (2.1), all the elementary normal forces 
point approximately toward the bearing center, so the bearing force components are 


b per Ep 
Fz =2 sing Rdg dz, Fy =2 cos g Rdg dz. (2.25) 
č die p sin p rag 7 rae ee p cos pọ KAP 


In anticipation of Chapter 3, the integrals in Eq. (2.25) assume the pressure field to be 
symmetrical in the axial direction within each bearing and consider only half the gap. 
Finally, it is emphasized that, in contrast to oil-lubricated bearings (see, e.g., [Sze10]), 
there is no need of neglecting subambient pressure in Eq. (2.25) since cavitation is 
not an issue with mostly gaseous lubricants. Similarly, it seems redundant to imitate 
foil lift-off by such a procedure in combination with a compliant foil model [LBS17b]. 
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3 Lubrication Theory for 
Phase-Changing Fluids 


The theory of fluid film lubrication, which goes back to British physicist and engineer 
Osborne Reynolds (1842-1912), allows for the calculation of pressure distributions in 
thin lubrication gaps on the basis of suitably simplified equations of fluid mechanics. 
Since some of the originally made assumptions must be questioned with regard to 
the non-ideal gas behavior in refrigerant-lubricated GFBs, this chapter introduces 
thermodynamic properties to enable the derivation of a modified Reynolds equation. 
Knowing that the behavior of refrigerants is quite sensitive to temperature variations 
in view of possible vapor-liquid phase transitions, an adapted energy equation for 
the prediction of temperature distributions is derived afterward in direct analogy. 
At the end of this chapter, having selected realistic boundary conditions for GFBs, 
all equations and constitutive laws can be combined into boundary value problems. 


3.1 Characterization of Volatile Refrigerants 


Before the actual fluid model can be properly derived, it is important to clarify some 
underlying assumptions and relationships regarding the thermodynamic behavior of 
compressible matter. Without claiming to be exhaustive, but in the most general way 
possible, the main concepts and ideas taken up later are therefore addressed below. 
Even though some of the considered correlations are illustrated by plots based on 
substance data of the hydrofluorocarbon refrigerant R-245fa, this does not mean that 
this section does not fully apply to other refrigerants or to similar fluids in general. 


3.1.1 Thermodynamic Fundamentals 


The fluid model is intended to describe the compressible flow of the lubricant within 
the narrowly confined space between the foil structure and the rotor shaft journal. 
Given that the permissible film thickness ho from Eq. (2.17) is often as small as 
a few micrometers, it appears essential to first assess the relevance of molecular 
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rarefaction effects before giving a continuous description of the fluid. When relating 
the length scale {1 of the available mean free path between the molecules to the 
critical macroscopic flow dimension, only a small characteristic Knudsen number 


Kn = fn <1 (3.1) 

ho 
actually guarantees the validity of the continuum hypothesis (see, e.g., [Kan+13]). 
Here, taking into account that the fluid compression in the narrowest zone of the 
lubrication gap typically reduces the mean free path to tens of nanometers [Pieoo], 
the expected Knudsen numbers are of the order of 107? and thus Eq. (3.1) is fulfilled. 


Under the assumption of a continuum, the current state of the moving fluid can be 
described by a number of continuous fields of intensive thermodynamic properties. 
This means that local values for pressure, density, temperature, or specific internal 
energy are obtained from averaging over a representative number of fluid molecules 
in the neighborhood of each point. Provided that the associated volume elements are 
sufficiently small to be treated as being homogeneous, it is consistent to additionally 
hypothesize that each of them is permanently in local thermodynamic equilibrium 
even though the fluid as a whole is not necessarily equilibrated (see, e.g., [Hueg8]). 
As a matter of fact, the general acceptance of local thermodynamic equilibrium even 
in highly dynamic systems is very closely related to why sufficiently small Knudsen 
numbers justify a fluid to be modeled as a continuum at all (see, e.g., [Kan+13]). 


When applying the laws of ordinary equilibrium thermodynamics, the state principle 
for compressible pure substances and non-reacting mixtures postulates that only 
two of the intensive properties can be chosen independently, while all remaining 
properties may be expressed through dependent state functions (see, e.g., [Mor+14)). 
In the pioneering work of Garcia [Garı2], pressure and temperature are identified as 
suitable primary variables for the analysis of refrigerant-lubricated GFBs in stationary 
operation. Especially when dealing with instationary compressible flow, however, 
it appears more natural to retain density instead of pressure as an unknown in the 
equations of fluid mechanics (see, e.g., [FPS20]). Consequently, the density field p 
and the temperature field Ÿ are selected hereafter as the two primary variables in 
terms of which the governing equations of the fluid model are to be formulated. 
As seen later in Sections 3.2 and 3.3, this choice is extremely advantageous because 
the respective time derivatives of density and temperature appear natively in the 
considered equations, which is not the case for the time derivative of pressure. 


For getting a qualitative impression of thermodynamic characteristics of a substance, 
the collectivity of possible equilibrium states is usually visualized as curved surfaces 
in three-dimensional phase diagrams that relate any three of the intensive properties. 
To obtain a p—p— phase diagram as shown in Fig. 3.1 for the refrigerant R-245fa, 
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Figure 3.1: Three-dimensional p—o-—0 phase diagram for the hydrofluorocarbon refrigerant R-245fa, 
which gives a qualitative impression of thermodynamic characteristics with possible phase transitions. 
According to the plotted data, the critical temperature is found to be approximately V. ~ 427 K (154°C). 


the dependent pressure p is plotted against the density pọ and the temperature 0 
by retrieving data from the open-source fluid property library CoolProp [Bel+14]. 
After fixing the thermodynamic state via density and temperature, such a diagram 
unambiguously indicates number and type of occurring phases (see, e.g., [Mor+14]). 
Here, according to the labels in the plot, it is distinguished between a vapor region, 
a liquid region, and a vapor-liquid region where both phases coexist in equilibrium. 
The so-called saturation dome is bounded by a saturated vapor line (dashed in gray) 
and by a saturated liquid line (dot-dashed in gray) that meet at some critical point. 


When following the exemplary isothermal lines (in magenta), it becomes obvious 
from the phase diagram in Fig. 3.1 that pressure increases monotonically with density, 
but never strictly monotonically below the critical temperature © of the refrigerant. 
Instead, all of the subcritical isotherms run horizontally through the saturation dome, 
where they remain on the level of particular isobaric lines (in green) that represent 
resulting saturation pressures at the respective temperatures (see, e.g., [Mor+14]). 
For the compressed refrigerant in a GFB, the practical consequence of a beginning 
phase transition somewhere in the lubrication gap is that the local pressure cannot 
increase further until either condensation is completed or the local temperature rises. 
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3.1.2 Formulation of Equations of State 


Among the abovementioned state functions that relate the two selected independent 
properties density and temperature to any other intensive property, some become 
particularly relevant in the following. On the one hand, a thermal equation of state 


p = pP, 0) (3.2) 


later allows for an elimination of the pressure p in the differential equations to be 
solved according to Section 3.2, but also for a subsequent calculation of the pressure 
distribution (see, e.g., [FPS20]). On the other hand, knowing that the fundamental 
equations of fluid mechanics describe transport mechanisms related to the specific 
internal energy field E rather than the evolution of temperature distributions, it is 
only a caloric equation of state 

E = E(p, 8) (3-3) 
that actually brings the temperature 0 into play while eliminating energy expressions. 


Despite the simple form of Eqs. (3.2) and (3.3), establishing such equations of state 
for real fluids is far from trivial and very often it is necessary to include numerous 
empirical parameters in models that satisfactorily reproduce experimental findings. 
As one of the most sophisticated descriptions of the considered refrigerant R-245fa, 
the recent approach by Akasaka et al. [AZL15] provides a very accurate model that 
proceeds via the Helmholtz free energy to deliver thermodynamic state functions. 
Based on the established equations of state, which are found sufficiently well-behaved 
throughout the entire range of interest, it is convenient to simply construct look-up 
tables that store a representative number of pre-calculated discrete function values. 
When evaluating data points in between by a suitable interpolation method [MHo1], 
a tabular formulation can be expected to compete with the modified Peng—Robinson 
equation of state proposed by Garcia [Gar12] in terms of both precision and efficiency. 


Using the tabular Helmholtz free energy formulation, which is available in CoolProp, 
the plot in Fig. 3.2a pictures Eq. (3.2) for R-245fa and shows pressure as a function of 
density and temperature along some isotherms (in magenta) and isochores (in blue). 
Compared to the phase diagram in Fig. 3.1, the considered density and temperature 
domains cover only such states that appear relevant for refrigerant-lubricated GFBs. 
Hence, the given representation focuses on the vapor region (V) and on the transition 
into the vapor-liquid region (V-L) across the saturated vapor line (dashed in gray), 
with a dot marking the state under normal conditions (i.e., p = 1 atm, 0 = 293.15 K). 
In a similar way, just rotated around the vertical axis for better visibility, the plot in 
Fig. 3.2b depicts the internal energy and therefore visualizes Eq. (3.3) for R-245fa. 
While the energy decreases only slightly with increasing density and with decreasing 
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temperature in the vapor region (V), it drops sharply immediately after entering the 
vapor-liquid region (V-L) due to a beginning release of latent heat of condensation. 


In the context of fluid mechanics, in which changes in internal energy are more of 
interest than absolute levels, it appears particularly useful to transfer the general state 
function from Eq. (3.3) to the corresponding exact differential (see, e.g., [Mor+14]) 


dE dE 
dE @ dd (&) do (3.4a) 
p v 
_ [9E p 1[ op 
(5) dd 02 h (3) | do (3.4b) 
p p 
= cy dd + ha — ay) do. (3.40) 


Having a closer look at the reformulated differential in Eq. (3.4b), which results from 
Eq. (3.4a) after a short calculation that uses the specific entropy (see, e.g., [Mor+14]), 
the thermal equation of state from Eq. (3.2) allows to exactly determine (0E/ 0p) 9. 
Since this form of the differential is extensively used in Section 3.3, an abbreviated 
notation is introduced in Eq. (3.4c), where the specific isochoric heat capacity field 


d 
Cv = Cv(p, v) = 5g Ee v) (3-5) 


is finally a function of the selected independent properties density and temperature, 
as is the isochoric compressibility coefficient field 


1 ð 
Ay = ay (p, 0) = EN 59 P(e, 8): (3.6) 


Later, this seemingly arbitrary linking of thermal and caloric equations of state leads 
to a particularly elegant form of the internal energy equation for compressible fluids. 


For the refrigerant R-245fa, the heat capacity is calculated by partially deriving* the 
internal energy that is obtained from the tabular Helmholtz free energy formulation. 
The plot in Fig. 3.2c visualizes Eq. (3.5) and reveals a discontinuity along the saturated 
vapor line (dashed in gray), which corresponds in Fig. 3.2b to the plotted sharp edge. 
In the vapor-liquid region (V-L), this implies from a physical point of view that cy 
acts as an apparent heat capacity that also considers the latent heat of condensation. 
Analogously, the plot in Fig. 3.2d pictures Eq. (3.6) for R-245fa based on the derivative 
of pressure, so that the sharp edge shown in Fig. 3.2a results in another discontinuity. 


t For practical reasons, the partial derivatives of thermal and caloric equations of state are approximated 
by applying a central finite-difference scheme to the tabular Helmholtz free energy formulation. 
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Figure 3.2: Plots of thermodynamic state functions for R-245fa from Helmholtz free energy formulation, 
with the saturated vapor line (dashed in gray) separating the vapor (V) from the vapor-liquid (V-L) region. 
For better orientation, a black dot marks the state under normal conditions p = 1 atm, 0 = 293.15 K (20°C). 
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Here, this expresses the fact that the volume drops substantially during condensation, 
which then causes a, as an apparent compressibility coefficient to grow abruptly. 
Later, some smoothing is applied to make state functions continuously differentiable. 


After all, there are so-called transport properties such as the dynamic viscosity field y, 
which are sometimes not considered thermodynamic properties in the strictest sense, 
but can still be represented as unique functions of the current state (see, e.g., [Obe6o)). 
Based on this finding, it is generally possible to assume a viscous equation of state 


u = (p, 0) (3.7) 


unless dealing with non-Newtonian fluids [Spigo], which is beyond the scope anyway. 
As Bell et al. [BL16] suggest for R-245fa and other refrigerants, there is a correlation 
between the residual entropy from a Helmholtz free energy formulation and viscosity. 
With regard to the representation of Eq. (3.7) that is given by the plot in Fig. 3.2e, 
the most important rise in viscosity occurs in the mixed vapor-liquid region (V-L). 
In view of the lever rule, this is due to a much higher viscosity of the liquid phase, 
which has a growing influence on the overall viscosity as condensation progresses. 


In thermodynamics, the term vapor quality usually refers to an intensive property 
that is defined only within the saturation dome and describes the mass percentage 
of vapor in a vapor-liquid mixture according to the lever rule (see, e.g., [Mor+14]). 
For subcritical states, this definition can be extended to a generalized vapor quality 


w = w(p,®) (3-8) 


that becomes w = 100% for a superheated vapor and w = 0% for a subcooled liquid. 
Even though vapor quality then loses its property as thermodynamic state function, 
this is unproblematic as long as density and temperature remain primary variables. 
On the whole, this means that Eq. (3.8) as plotted in Fig. 3.2f follows directly from 
the previously determined shape of the saturation dome for the refrigerant R-245fa. 
Later on, this property is not explicitly included in any of the governing equations, 
but rather serves as a measure to characterize the phase condition of resulting states. 


3.1.3 Validity of Perfect Gas Assumption 


A closer look at the plots in Fig. 3.2 quickly reveals that deep in the vapor region (V) 
and thus far from the saturated vapor line (dashed in gray), all thermodynamic state 
functions appear linear or even constant with respect to density and temperature. 
Hence, the simple equations of state that emerge from the idea of a perfect, ideal gas 
seem justified there and can be used later for the purpose of model comparisons, 
knowing that these are exactly the equations commonly used for air-lubricated GFBs. 
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Table 3.1: Constants in the perfect gas model of R-245fa. 


Description Symbol Value Unit 
Specific gas constant Rs 620 J/(kg-K) 
Specific heat capacity c3 803 J/(kg-K) 
Heat of condensation AE° 177 J/g 
Viscosity w 11.6 yPa-s 


Assuming the pressure to be proportional to the product of density and temperature 
with the specific gas constant Rs as the proportionality constant (see, e.g., [Mor+14]), 
the thermal equation of state from Eq. (3.2) becomes the well-known ideal gas law 


p = p(p, 0) = pRsd. (3.9) 


Inserting Eq. (3.9) and its partial derivative with respect to temperature into the 
definition given by Eq. (3.6), the compressibility coefficient takes the compact form 


1 

3: (3.10) 
Next, the internal energy of a perfect gas being linear with respect to temperature 
but independent of density, the caloric equation of state from Eq. (3.3) is written as 


E=E(8) =cj0+AE°. (3.11) 


In Eq. (3.11), the constant slope is directly identified with the specific heat capacity 

cy = 08 (3.12) 
stated by Eq. (3.5), while the constant AE° represents the latent heat of condensation. 
To keep the level of complexity consistent, it is assumed hereafter that the viscosity 


n=w (3.13) 
and therefore the viscous equation of state from Eq. (3.7) are also constant properties. 
Except for Rs, which is an actual constant for a given substance, all other constants 
in Table 3.1 are calibrated to match values from CoolProp under normal conditions. 
Lastly, with regard to the generalized vapor quality from Eq. (3.8), setting w = 100 % 
is a trivial statement since there is no liquid phase by definition of a perfect gas. 


To demonstrate the limits of the perfect gas model, the two plots in Fig. 3.3 illustrate 
the relative errors that are made in relation to the Helmholtz free energy formulation. 
In the vapor region (V), the perfect gas reproduces all the properties quite precisely, 
not only in Fig. 3.3a near the reference point at the normal temperature 0 = 293.15 K, 
but also in Fig. 3.3b at the temperature 0 = 373.15 K with still acceptable deviations. 
In the vapor-liquid region (V-L) in Fig. 3.3a, however, the clearly failing perfect gas 
shows the need for realistic equations of state to model refrigerant-lubricated GFBs. 
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Figure 3.3: Plots of the relative errors in various properties that are made under the perfect gas assumption 
in relation to the Helmholtz free energy formulation. For the viscosity, e.g., this is |[u (p, 8) — u°]/u(p, ®)|. 


3-2 Derivation of Modified Reynolds Equation 


Regarding the fluid as continuous matter on the basis of the small Knudsen number 
from Eq. (3.1), the flow is described by the usual equations of fluid mechanics, namely 
the Navier-Stokes equations, the continuity equation, and the energy equation. 
Since these coupled equations can neither be solved analytically nor be evaluated 
numerically without considerable computational effort, it is often helpful to simplify 
them as much as possible for a specific application. For this reason, the thin-film 
assumption from Eq. (2.1) is used below and in Section 3.3 to ultimately derive a 
lubricant model that is rather simple but still sufficiently accurate for most purposes. 


3.2.1 Approximated Navier-Stokes Equations 


First of all, the Navier-Stokes equations for compressible Newtonian fluids are 
written in the general vector form 


radu + grad'u 


poe = pf — grad p + grad (À div u) + div | 


which can be found in a majority of textbooks on fluid mechanics (see, e.g., [Hue98]). 
The differential operator 


D d 

De (5 gn (3.15) 
here denotes the material derivative with local and convective time rates of change for 
a physical field in Eulerian description. The given partial differential equation (PDE) 
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serves as the most common mathematical model for relating the pressure field p 
and the density field p to the flow velocity vector field u in a moving isotropic fluid. 
More specifically, the underlying principle of conservation of momentum states that 
the inertial forces on the left-hand side of Eq. (3.14) are balanced by the sum of body 
forces, pressure forces, and viscous forces. In the first term on the right-hand side, 
the body force (per unit mass) vector field f needs to account only for gravitational 
acceleration because the fluid is observed from an inertial frame of reference, so that 
no additional centrifugal or Coriolis accelerations appear. Moreover, as expressed by 
the last two terms on the right-hand side, viscous stresses in a Newtonian fluid are 
linearly proportional to strain rates and may depend on the dynamic viscosity field y 
as well as on the volume viscosity field A in case of compressible behavior. Although 
the considered GFBs do not operate with ambient air, but still with refrigerants of 
comparable rheology, non-Newtonian effects as only reported for some rarely used 
multi-component lubricants [DE83] are expected to play an entirely subordinate role. 


When expressing Eq. (3.14) in terms of the Cartesian coordinates (x,y,z) formerly 
introduced for the lubrication gap, the existence of two greatly different length scales 
justifies considerable simplifications of the involved differential operators. With the 
ratio between these characteristic length scales resulting in the small parameter € 
according to Eq. (2.1), it is commonly argued in fluid film lubrication theory that 
ratios of the same order of magnitude must hold between the respective flow velocity 
components (see, e.g., [Sze1o]). Having said this, all the coordinates (including time t) 
as well as the flow velocity components ux, Uy, and uz can be rescaled to normalized 

variables of the order of unity (denoted hereafter with a tilde) by substituting 
x = lož, y = Eloÿ, z = loz, t= 2 ; (3.16a) 

wo 

Ux = lowollx, Uy = Elowofy, Uz = lowoüz. (3.16b) 


Based on the angular frequency wo from Eq. (2.11), the characteristic time scale wg 1 
in Eq. (3.16a) and the characteristic velocity scale {owo in Eq. (3.16b) describe a flow 
dominated by the tangential boundary motion of the rotating shaft journal [Bau18]. 
In an analogous manner, the dynamic viscosity y, the volume viscosity A, the fluid 
density p, the fluid pressure p, as well as the body force components fx, fy, and fz 
are transformed to normalized variables by setting 


g x sy wo. 
u = Hol, À = pod, p = pop, p= F = 5, (3-17a) 


fx = gfe fy = shy fz = gf- (3.17b) 


The characteristic viscosity ig is representative for both viscosity coefficients u and A 
in Eq. (3.17a) since values of the same order of magnitude can be expected in the 
light of Stokes’ hypothesis 2u + 3A = 0 for monatomic gases, which holds at least 
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approximately in the considered case of polyatomic gases (see, e.g., [Hueg8]). Besides 
the introduction of a characteristic density pọ in Eq. (3.174), it appears particularly 
noteworthy that a viscous pressure scale jigwoe * instead of an inviscid pressure 
scale has to be adopted in order to retain pressure as an unknown field quantity in the 
reduced PDE system, which would be overdetermined otherwise (see, e.g., [Sze10]). 
Finally, under the aforementioned assumption of body forces resulting only from 
gravity in the considered problem, the constant of gravitational acceleration g in 
Eq. (3.17b) gives the correct order of magnitude for all three body force components. 


Conducting an order-of-magnitude estimation for the normalized Navier-Stokes 
equations in Cartesian coordinates, most of the viscous stress terms are of orders 
greater or equal to £? and thus appear negligible compared to unity in virtue of the 
thin-film assumption from Eq. (2.1), which yields the reduced PDE system 


dil dil oû dil &Re_- op af _oû 
2 re x 3 x 4 x fe on NENS > P - OUx 
. Rea D Or og op | RP 3g + 3g (r ) 6.188) 


g oy 
dü oû dü ou Re ~ op 
An» = Yi ee Sy CY. PR E p 
€ Rea y tx À hy 37 iz | => Pfy ag’ (3.18b) 
dü oû dü oû &Re_- op Of _dii 
2 x z i z 5 z = ZN 2. ~ P ~ OUz 
€ Rea JE Üx SE tly ET Uz a5 | E Õfz a5 * 3 (r a ) (3.18c) 


Among the remaining terms, the relative importance of inertial forces and of body 
forces is then determined by a characteristic Reynolds number and a characteristic 
Froude number that are, respectively, defined as 


2 | 
Re = nun j Fr = wọ í A (3.19) 


Looking at typical values for the involved physical quantities such as introduced in 
Chapter 6, the scaled Reynolds number & Re in both Eqs. (3.18a) and (3.18c) ranges 
between 10~2 and 10~!, whereas ¢2 Re Fr? is located in the region of 10° to 1074. 
With the exception of excessively high rotational speeds, this confirms the essential 
hypothesis by Reynolds [Rey86] that inertial forces and body forces acting upon a 
thin fluid film appear insignificant compared to pressure forces and viscous forces. 


Finally, the coefficients e Re and £3 Re Fr in Eq. (3.18b) being smaller than even 2, 


the back-transformed expression that remains from the original PDE simply writes 
ð 
Sri 
dy 


and suggests that the pressure field p is independent of the cross-film coordinate y. 
In classical lubrication theory, it is additionally hypothesized that not only pressure 


(3.20) 
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but also density, viscosity, temperature, and other properties are constant across a 
sufficiently thin fluid film (see, e.g., [Cop49]). This simplification is overcome, at least 
in theory, with the generalized Reynolds equation proposed by Dowson [Dow62], 
which accounts for possible cross-film variations of density and viscosity. In practice, 
however, it can be extremely cumbersome or even impossible to rigorously describe 
the transient behavior of an entire bearing-rotor system starting from such a partial 
integro-differential equation. In the literature, cross-film variations of properties are 
considered almost exclusively in very detailed steady-state analyses [Gar12, MLS16], 
whereas comparable rotordynamic studies [Bou17, LMS18] understand especially 
density as cross-film averaged quantity. As recently outlined by Lehn et al. [LMS18], 
the generalized Reynolds equation has of course its justification, but rather when 
applied to oil-lubricated bearings because the magnitude of non-uniformities across 
gas films is often insignificant. With this in mind, the present thesis supposes that 
the density field is a priori independent of the cross-film coordinate in analogy to 
what Eq. (3.20) indicates for pressure. Having said this, the thermodynamic state 
principle mentioned in Section 3.1.1 requires that, since the two properties pressure 
and density do not depend on y, any other properties such as viscosity do neither. 
Setting accordingly 


op _ Ou _ 
T =0, où = 0, (3.21) 
the back-transformed terms remaining from Eqs. (3.18a) and (3.18c) simplify to give 
op a Oux\ _ ur Gon) 
dx doy K dy K oy? ” ue 
op a duz \ _ u; 
ae (r dy ) u TA (3.22b) 


To actually obtain flow velocity profiles depending on prevailing pressure gradients, 
both Eggs. (3.22a) and (3.22b) must be integrated twice with respect to y and can then 
be adapted to boundary conditions that are related to the motion of adjacent walls. 


On the inner fluid film boundary, which is found at y = —h, according to Eq. (2.14), 
the tangential component of the wall velocity is dominated by the rotation rather than 
by the translation of the shaft journal and points approximately in the e,-direction 
for a sufficiently thin film (see, e.g., [Sze10]). Together with the simplifications that 
any tangential motion of the foil structure on the outer boundary at y = —q appears 
negligible [Bau18] and that neither of the bounding walls moves axially, this yields 


à = Rwo, Ux 
y=—he 


= (3.232) 
ya-q 


_ =0. (3.23b) 
Yad 


Ux 


= 0, Uz 
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Possible refinements of such no-slip boundary conditions, which could be extended 
to consider slip flow as proposed by Burgdorfer [Bur59], are beyond the scope of 
the present thesis. Nevertheless, it should be kept in mind that there are possible 
applications with elevated Knudsen numbers, in which no-slip boundary conditions 
might become inaccurate because of beginning molecular rarefaction effects [HD83]. 
Here, keeping simply with Eq. (3.23a), the circumferential flow velocity field 


19 + 
Ux = > u (hs y) (q y) DE +) (3.24) 


matches the prescribed velocity Rwo on the inner fluid film boundary at y = —h, 
and vanishes on the outer boundary at y = —q, whereas the axial flow velocity field 


19 
as, = (ha + y)(q+y) (3.25) 


satisfies Eq. (3.23b) and is zero at both boundaries. Generally speaking, the effective 
circumferential flow that is described by Eq. (3.24) results from the superposition of 
a pressure-driven Poiseuille flow and a shear-driven Couette flow, whereas the axial 
flow described by Eq. (3.25) corresponds to a purely pressure-driven Poiseuille flow. 


3.2.2 Integration of Continuity Equation 


To proceed with the derivation of an appropriate lubricant model, the Navier-Stokes 
equations from Eq. (3.14) are completed by the continuity equation for compressible 
fluids in the general form 


Dp naa 
Dr + p divu =0, (3.26) 


which again is found in all common textbooks on fluid mechanics (see, e.g., [Hue98]). 


With the definition of the material derivative from Eq. (3.15) and with a generalized 
product rule, the PDE from Eq. (3.26) can be stated in Cartesian coordinates as 


ar + an (Pra) + (ou) + ge (Pre) =0 627 


where all terms are of the same order of magnitude if the characteristic length and 
velocity scales from Eqs. (3.16a) and (3.16b) hold. Since no further information about 
the radial flow velocity uy can here be obtained from the Navier-Stokes equations, 
an essential step consists in the integration of Eq. (3.27) across the fluid film, giving 


Ft fy gale) du + (om) 


FA 


-1 à 
EP + L 32 (puz) dy = 0. (3.28) 
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In Eq. (3.28), the unknown u, appears only at the boundaries of the fluid film, where 
it can be directly related to the motion of the shaft journal and of the foil structure. 


The normal components of the wall velocities at the respective boundaries, which 
point both approximately in the ey-direction for a sufficiently thin film, are given by 


dhs ah aq 


= R = SS : 
“y y=-h; ot a ax ” "y y=-q ot (3-29) 


and consider the rotation and the translation of the shaft journal (see, e.g., [Sze1o]) 
as well as the deformation rate of the foil structure [Bau18]. Since the density field p 
is supposed to be constant across the film according to Eq. (3.21), the third term on 
the left-hand side of Eq. (3.28) can directly be expressed with the help of Eq. (3.29) as 


(em) 


í oq Ihr Ons 
-hx = e(l 7 uy n) Pa, tea, PRV 630) 


With Leibniz’ rule for differentiation under the integral sign, which is obviously 
applicable to each of the three integrals on the left-hand side of Eq. (3.28), it comes 


-139 al m ə əh, 
Wale) Me. A 

-41 à 0 —4 oh 
Late) dy = AG a va) a GBT) 


-1 à d -4 
I È = (ouz) dy = > (e a i say), (3310) 
where density is constant across the film in accordance with Eq. (3.21) and velocities 
at the integration limits are obtained from Eqs. (3.23a) and (3.23b), respectively. 


When finally rewriting Eq. (3.28) using these conversions, all terms from Eq. (3.30) 
cancel out with the corresponding terms in Eqs. (3.31a) and (3.31b), so it only remains 


ð -4 d ð —4 d 0 4 Avl=o 

ale], ©) +35 p |, ay ae p f way) =0. (3.32) 
Remembering the relation h = h, — q from Eq. (2.13c), the trivial integral in the first 
term of Eq. (3.32) simply yields the fluid film thickness field h itself. Additionally, the 


two remaining integrals can be compactly expressed as hü, and hū, after calculating 
from Eggs. (3.24) and (3.25) the cross-film averaged bulk flow velocity fields 


EEE de h? op Roo _ - i 

G a MOY ig a Pe (3-33a) 
1 h2 ap 

=, | HU nn (3.33b) 
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According to Eq. (3.33a), the expression for the circumferential velocity field can be 
split into a first term üy p representing the Poiseuille flow and a second term ily c 
representing the Couette flow. In the case of a pure Poiseuille flow as in Eq. (3.33b) 
for the axial direction, such a decomposition of the velocity field is of course needless. 


To obtain a more descriptive and more common form of the resulting shortform PDE 
d ð 2 d Ni 
3 (ot) + È (as) + (m) =0 gan 
the bulk flow velocities from Eqs. (3.33a) and (3.33b) are explicitly inserted to yield 
d Rwo 9 Ə / ph? op d { ph? ap 
h) 4 h) = 5 
ot (e ) 2 ox (e ) ox (E ox dz \ 124 dz (3:35) 
or, making use of the arc-length transformation x = Rọ from Eq. (2.6), to finally give 
ə wg ð _ 10/phop\ , a { ph° op 
ot (er) 20% (er) R209 E og) ` əz \ 12u əz J` (3.36) 


Based on the fundamental work of Reynolds [Rey86] about incompressible lubricants, 
this PDE is subsequently referred to as the Reynolds equation for compressible fluids. 


3.2.3 Closure of Coupled Equation System 


In the most general form of Eq. (3.36) to be considered, the involved field quantities 


p = p(p,zt) p = p(9,2,8), u = p(p,2,t), h=h(g,t) (3-37) 


may depend on the time t as well as on the circumferential coordinate p and, with the 
exception of the fluid film thickness in the last expression, on the axial coordinate z. 
With regard to the number of unknown fields in Eq. (3.37), the Reynolds equation 
needs to be closed by some constitutive laws for the lubricant, which are a thermal 
equation of state to eliminate the pressure field p(@,z,t) according to Eq. (3.2) and 
a viscous equation of state as in Eq. (3.7) to express the viscosity field u(9,2z,t). 
Concerning the refrigerant R-245fa, for which the perfect gas assumption according 
to Eqs. (3.9) and (3.13) is clearly insufficient, the consideration of equations of state 
as shown in Figs. 3.2a and 3.2e adds an important nonlinearity to the resulting PDE. 
Moreover, as explained in Section 2.3, it is finally the insertion of h(g,t) from either 
Eq. (2.13a) or Eq. (2.13b) into Eq. (3.36) that introduces the coupling quantities from 
the foil structure submodel and from the rotor submodel into the fluid submodel. 
Provided the coupling quantities are adequately described by the equations of the 
other two submodels, the Reynolds equation is ultimately formulated only in terms 
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of the selected primary variables, which are the density field o(g,z,t) and the 
temperature field 0(g,z,t). Of course, a PDE that physically describes the evolution 
of the latter is missing up to this point, which is the motivation for additionally 
deriving a temperature equation in Section 3.3 based on the internal energy. 


In connection with the Reynolds number? from Eq. (3.19) and notwithstanding 
the assumption of negligible inertial forces, it remains to clarify the role of fluid 
instability mechanisms that lead to the formation of Taylor vortices or turbulent flow. 
Already in the pioneering theoretical investigation of turbulent gas lubrication by 
Constantinescu [Con64], it actually turns out that the influence of such phenomena is 
only of a quantitative nature and less significant in gas bearings than in oil bearings. 
For this reason and for the sake of simplicity, the widely used class of zero-equation 
turbulence models leaves governing PDEs such as the Reynolds equation basically 
unchanged? and introduces an additional turbulent eddy viscosity instead, which is 
algebraically related to quantities such as flow velocity gradients (see, e.g., [FPS20]). 
With such an approach, more precisely with the Ng-Pan-Elrod model [NP65, EN67], 
the numerical investigation by Garcia [Gar12] reveals that turbulence does occur in 
refrigerant-lubricated GFBs, but typically at locations that are largely insignificant 
for their overall performance. Although effects on bearing-rotor systems are thus 
limited, it is noted that the Reynolds equation in Eq. (3.36) principally proves suitable 
for a qualitative, but also for a quantitative description of turbulent flow provided 
the viscosity field is altered appropriately. Unsurprisingly, such a model extension 
would be relatively straightforward [Gar12, Bou17], but is not pursued further here 
to keep the number of empirical parameters small when analyzing qualitative trends. 


3-3 Derivation of Thin-Film Temperature Equation 


While it might be sufficient to generally assume isothermal conditions for the analysis 
of GFBs lubricated with ambient air, this is certainly not a reasonable supposition 
when dealing with refrigerant-lubricated GFBs [Garı2]. Consequently, the continuum 
hypothesis according to Eq. (3.1) and the thin-film assumption from Eq. (2.1) are 
used hereafter to derive a temperature equation similar to the Reynolds equation. 


? To be precise, the characteristic Reynolds number must be completed by a characteristic Taylor number 
to correctly distinguish between parallel flow instabilities and centrifugal instabilities (see, e.g., [Sze1o]). 


3 The form of the Reynolds equation for turbulent flow is unchanged since the underlying Reynolds- 
averaged Navier-Stokes equations also have the same form as the original equations (see, e.g., [FPS20]). 


34 


3.3 Derivation of Thin-Film Temperature Equation 


3.3.1 Approximated Internal Energy Equation 


According to usual textbooks on fluid mechanics (see, e.g., [Hueg8]), the combination 
of the Navier-Stokes equations from Eq. (3.14) with the continuity equation from 
Eq. (3.26) is typically accompanied by the internal energy equation for compressible 
Newtonian fluids in the general form 

2. 

DE 


DE gradu + grad'u 
P Di 


2 


= —pdivu — div p + A(divu)* +24 (3.38) 


Using the material derivative from Eq. (3.15), this PDE expresses the principle of 
conservation of energy by relating changes within the specific internal energy field E 
to different energy sources and energy sinks. While the first term of the right-hand 
side of Eq. (3.38) is associated with compression work, the second term contains the 
heat flux vector field @ and is therefore supposed to represent thermal conduction. 
Most importantly, the last two terms* correspond to viscous dissipation work, which 
depends on both the dynamic viscosity field u and the volume viscosity field A. 


Transferring the reformulated exact differential dE that is given by Eq. (3.4c) into the 
corresponding material derivative DE/ Dt and then substituting Dp/ Dt = -pdivu 
by virtue of the continuity equation from Eq. (3.26), it comes immediately 


DE Do p Do 
p Dt pev Dt l p (1 ayt) Dt (3:39) 
Dé : 
=p p div u + ay 0p div u. (3.39b) 


With this relationship, the internal energy equation from Eq. (3.38) takes a convenient 

temperature-explicit form to describe changes within the temperature field 0 by 

Tr |\2 

gradu + grad u 
2 


Do 
pcv ir = bp div u — div p + A(div u) +24 (3:40) 


This reformulated PDE now involves the specific isochoric heat capacity field cy 
on the left-hand side and the isochoric compressibility coefficient field ay in the 
first term on the right-hand side, whereas contributions due to thermal conduction 
and dissipation work remain unchanged. An alternative temperature-explicit form 
of this PDE is actually shown in [LBS19], where the specific isobaric heat capacity 
and the isobaric compressibility coefficient occur instead of isochoric properties. 


4 The mathematical expression for the viscous dissipation work makes use of the Frobenius norm ||.. .||, 
which is calculated as the square root of the sum of the squared tensor components (see, e.g., [GV13]). 
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However, such a derivation via the specific enthalpy leads to an appearance of the 
time derivative of pressure, which cannot be substituted by virtue of the continuity 
equation and thus makes the solution of transient problems much less intuitive. 


In the following step, in which the PDE from Eq. (3.40) is finally represented in the 
Cartesian coordinates (x,y,z) that describe the lubrication gap, the substitutions 
from Eggs. (3.16a) and (3.16b) allow again for a normalization of all coordinates and 
flow velocity components. Moreover, the various scalar fields of fluid properties as 
well as the heat flux components #;, y, and ¢z can be transformed to normalized 
variables (denoted again with a tilde) by means of Eq. (3.17a) and by setting 


2 
= Howo à E“ PoCo ~ 
= ® = v = . 
Cy = Coty, ieee ay RR dy, (3-41a) 
av Avo ~ u Ado ~ er Avo ~ Ga 
px = Ho, Par $y = COMO Pr fz = copo by Qz. 3-4 


After introducing in Eq. (3.41a) a characteristic specific heat capacity co, it is natural to 
construct the temperature scale uowoe 7p ‘co! and, reciprocally, the compressibility 
coefficient scale, in analogy to the adopted viscous pressure scale (see, e.g., [Cop49]). 
In Eq. (3.41b), assuming that heat fluxes generally obey Fourier’s law p « — grad 0, 
the constant of proportionality is of the order of couo according to the kinetic theory 
of gases (see, e.g., [Ken38]) and the temperature gradient may be assessed by relating 
a characteristic temperature difference A to the respective length scales {5 and ef. 


From a subsequent order-of-magnitude estimation for the normalized internal energy 
equation in Cartesian coordinates, it turns out that most of the terms related to 
viscous dissipation work are of orders greater or equal to £? and can thus be neglected. 
In doing so, it remains only the reduced PDE 


ee (3 3 ab | 2) 3 de dily =) 
ply + üx + ily + Üz = Add — + — + 


of 
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in which the relative importance of thermal conduction can directly be judged by a 
characteristic Eckert number defined as 


Gui 
an (3-43) 


Since Ec”! is expected of the order of unity when inserting typical parameter values 
according to Chapter 6, the cross-film heat flux must clearly be retained in Eq. (3.42), 
while heat fluxes along the film are smaller by the factor of £? and appear negligible. 
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On account of the state principle and with an intent to describe internal energy at the 
same level of model complexity as pressure and density in Eqs. (3.20) and (3.21), it is 
assumed hereafter that temperature, specific heat capacity, and the compressibility 
coefficient are all independent of the cross-film coordinate y, in other words 


dP _ 4 dv _ 4 y 


Even though this simplification seems entirely justified for investigating qualitative 
effects that concern the overall system, it should of course not be ignored that the 
accuracy of quantitative predictions can still be improved by allowing temperature 
and other properties to vary across the fluid film [Gar12, Bou17]. In order not to 
overcomplicate the fluid model without necessity, however, the present thesis sticks 
to Eq. (3.44), especially since there is both analytical [Sny65] and numerical [Leh17] 
evidence that non-uniformities across gas films are much less significant than they 
are across liquid films. Ultimately, this makes the cross-film temperature gradient 
disappear on the left-hand side of Eq. (3.42) and the back-transformed PDE becomes 


[N {aux , Ay de 
Cpe kag az Pl ox ay | à 


The flow velocity components and the retained heat flux component may still depend 
on y as they are not intensive thermodynamic properties obeying the state principle. 


3.3.2 Integration of Internal Energy Equation 


As no further information is retrieved from any equation about the fields u, and ¢y, 
the most convenient remedy is to integrate Eq. (3.45) across the film, which leads to 


ue OO 1,00 1,90 d 
[es dt əx “oz Y 


za Oux 
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—4 duz 
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| -4 Ouy q du; 
TREFN x an fafa) dy (3.46) 


in the very same manner as previously seen in Eq. (3.28) for the continuity equation. 
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After integration, the unknown cross-film heat flux ¢, appears only at the boundaries 
of the fluid film, where it is prescribed by simple heat transfer boundary conditions 


ty], =—Bor(?- 0x), y| __, = Bos(®- 8%). (63-47) 


y=—hs 
Here, imagining the entire bearing-rotor system as a network of thermal resistances, 
heat fluxes from or into the environment are assumed proportional to the difference 
between the local fluid temperature 0 and the constant ambient temperature 9%. 
With this empirical approach (see, e.g., [CG15]), all contributions due to convection, 
conduction, and radiation are combined into two equivalent heat transfer coefficients, 
of which Bo refers to heat fluxes through the rotor shaft and Bos to heat fluxes 
through the foil structure and bearing sleeve. Even though this concept obviously 
represents a strong idealization of actual heat transfer mechanisms, the in-depth 
analysis of GFB temperatures by San Andrés et al. [SK10] shows that the complexity 
of such a model is largely sufficient to reproduce experimental results. Consequently, 
the actual difficulty arises rather from the correct determination of empirical heat 
transfer coefficients, which is of course not the primary interest of the present thesis. 
Finally, the first term in the third line of Eq. (3.46) can be compactly expressed as 


|, = 7 = di = Bo (0 = dx) (3:48) 


I y--h 


when summing up the heat transfer coefficients in Eq. (3.47) to give By = Bos + Bor- 


The further evaluation of the integrated energy equation turns out to be relatively 
simple because, owing to Eqs. (3.20), (3.21), and (3.44), all remaining field quantities 
apart from the components of the flow velocity are independent of the y-coordinate. 
Having this in mind and starting with the most obvious part, the integral describing 
temperature changes in the first line of Eq. (3.46) immediately becomes 


où [9 où [9 où [9 
(ya [mar (war) (3.49) 


Furthermore, the terms related to compression work in the second line of Eq. (3.46) 
possess a structure that is somewhat similar to the integrated continuity equation, 
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from which one benefits by slightly modifying Eqs. (3.30), (3.31b), and (3.31c) to get 


=q dux 0 —4 oh, 
Te ay Op Jx dy ay Op = 5 ux dy — aydpRwo Jx” (3.50a) 
= oq oh, oh, 
(avopuy) o = —ay0p Pr ay Op = 1 avðpRwo a, (3.50b) 
-q duz _ o [1 
= ay Op a dy ay Op z Ja uz dy. (3.50c) 


On closer inspection of the found expressions, all of Eqs. (3.49) and (3.50a) to (3.50c) 
can be stated in terms of the effective fluid film thickness h = h, — q from Eq. (2.13c) 
and of the bulk flow velocities ñy and ü; from Eqs. (3.33a) and (3.33b), respectively. 


Coming to the integrals of viscous dissipation work in the third line of Eq. (3.46), 
however, no such simplification seems possible without prior knowledge about the 
involved flow velocity fields. Consequently, it is only after calculating the squared 
derivatives of Eqs. (3.24) and (3.25) that carrying out the actual integration yields 
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Surprisingly, each summand in Eqs. (3.51a) and (3.51b) resembles the square of a 
particular bulk flow velocity contribution when setting h = h, — q and comparing 
the resulting expressions to Tp and Tc from Eq. (3.33a) and to #2 from Eq. (3.33b). 
Since there is no general relationship between the integrated square of the derivative 
of an arbitrary function and the squared integrals of its individual summands, this is 
only a neat side effect for flow velocity profiles given by second degree polynomials. 


Altogether, having substituted Eqs. (3.48), (3.49), (3.50a) to (3.50c), (3.51a), and (3.51b) 
into Eq. (3.46), the final temperature equation is given by the handy shortform PDE 
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Alternatively, when inserting the bulk flow velocities from Eqs. (3.33a) and (3.33b) 
and additionally exploiting the arc-length transformation x = Rp stated by Eq. (2.6), 
the temperature equation in Eq. (3.52) can be brought to the more descriptive form 
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Together with the particular form of the Reynolds equation in Eq. (3.36), the derived 
temperature equation governs the transient behavior of any compressible lubricant. 


3.3.3 Closure of Coupled Equation System 


In addition to the field quantities listed in Eq. (3.37), the PDE in Eq. (3.53) involves 


0 = 0(9,z,t), Cy = Cy(9,zZ,t), Ay = ty(@,zZ,t), (3-54) 


all of which depend once more on the time t as well as on the coordinates @ and z. 
Here, to close the temperature equation, exactly the same constitutive laws as for 
the closure of the Reynolds equation are required, which are a thermal equation of 
state and a viscous equation of state according to Eqs. (3.2) and (3.7), respectively. 
While the former already allows for an elimination of the compressibility coefficient 
field a, (g,z,t) owing to Eq. (3.6), an additional caloric equation of state as in Eq. (3.3) 
must be known to express the heat capacity field cy(ọ,z,t) as Eq. (3.5) stipulates. 
Moreover, the temperature equation is not only coupled to the Reynolds equation, 
but also to the other submodels via the insertion of h(ọ, t), which makes the included 
coupling quantities the only unknowns besides the fields p(g,z,t) and 8(,z,t). 


3.4 Statement of Boundary Value Problems 


Remembering that the Reynolds equation in Eq. (3.36) and the temperature equation 
in Eq. (3.53) are based on the thin-film assumption from Eq. (2.1), neither of them 
remains valid beyond the axial bearing edges or in the vicinity of foil fixation gaps. 
The PDEs are thus solved only within certain open and bounded domains N, C R?, 
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whereas additional constraints are provided on the corresponding boundaries 0Qn. 
Altogether, this creates one or more boundary value problems (BVPs) with solutions 
that belong to spaces of sufficiently many times continuously differentiable functions. 
At this point, it can already be anticipated that the considered number of pads Np 
determines not only the geometric shapes of (2, and 9(2,, but also the total number 
of individual BVPs into which the overall fluid problem is subsequently divided. 


3.4.1 Classification of Partial Differential Equations 


In the theory on linear PDEs of at most second order, the unambiguous classification 
into equations of either elliptic, parabolic, or hyperbolic type decides on appropriate 
boundary and initial conditions that lead to a well-posed problem (see, e.g., [Are+18]). 
For the nonlinear PDEs of fluid mechanics, on the other hand, such a distinction is 
much more delicate because the expected characteristics for these equations are often 
the result of combined elliptic, parabolic, and hyperbolic behavior (see, e.g., [FPS20]). 


With regard to the Reynolds equation in Eq. (3.36), it can in particular be observed 
that the convective terms on the left-hand side are hyperbolic in themselves, whereas 
the equally important viscous terms on the right-hand side have a parabolic character. 
The temperature equation in Eq. (3.53) is obviously dominated by convective terms 
of hyperbolic character, which means that only the few terms that are similar to the 
right-hand side of the Reynolds equation can be associated with parabolic behavior. 
In possibly existing regions of a flow where pressure gradients physically disappear, 
for example due to an isothermal condensation at prescribed vapor pressure [GS19], 
the parabolic contributions in both of the PDEs may locally even be completely lost. 


Besides differences in the choice of suitable solution methods, the crucial point is that 
well-posed parabolic problems require boundary conditions at all edges, whereas for 
hyperbolic problems it reveals unphysical to impose restrictions at outflow edges. 
To enable such a distinction, each boundary is formally partitioned according to 


90, 80, U INQn,o OO Rs (3-55) 


into three pairwise disjoint and not necessarily connected regions, namely an inflow 
boundary d02,:, an outflow boundary 0Qzy,o, and finally a symmetry boundary 902,,s- 
In the following, the assignment of edge points to either the inflow or the outflow 
boundary is made on the basis of the respective bulk flow velocity components, 
which implies that the partitioning of the boundaries is a priori unknown. 
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3.4.2 Identification of Domains and Boundaries 


In the case of a full-arc bearing with Np = 0, there is one single domain denoted by 
1 L 


which corresponds topologically to the lateral surface of a cylinder with height L/2. 
With the definition of the circumferential coordinate on the circle T! = R/(27Z), 
only 27t-periodic solutions are permitted as the schematic sketch in Fig. 3.4a indicates. 
In the axial direction, it is sufficient to consider only half the bearing length owing to 
the symmetry that results from identical conditions on both sides of the bearing and 
from supposing the film thickness in Eq. (2.13) to be independent of the z-coordinate. 
Accordingly, the symmetry boundary of a zero-pad GFB writes 


95 = T! x {0} (3.57) 


and is directly identified with the circumferential centerline of the bearing at z = 0. 
At the open bearing edge at z = L/2, the fluid has the possibility to either leave 
or enter the lubrication gap, which can be judged by whether the axial pressure 
gradient dp /dz makes the velocity component üz from Eq. (3.33b) positive or negative. 
Relying on this distinctive criterion, the mathematical expressions 


IM; = foz € T! x {5} 


IN, = foz ET! x 15} 


describe, respectively, the inflow and outflow boundary of a full-arc bearing without 
of course any possibility of fluid flowing through the non-existent foil fixation gaps. 


op 
a > ; 
De o}, (3.58a) 


op 
3z < o) (3.58b) 


Coming to the far more important case of a multi-pad bearing with Np > 1 pads, 
the situation is complicated by the fact that there are now Np unconnected domains 


L 
Q, = (Xn, Xn+4x) x (o, a n=1,...,Np, (3.59) 


which are separated by the Np foil fixation gaps as the sketch in Fig. 3.4b shows. 
Based on Eqs. (2.2) and (2.3), the circumferential coordinate g for the n-th domain 
must range within the leading gap position x, and the trailing gap position? Xn+4Ax. 
Given that the abstracted fixation gaps themselves have virtually no spatial extension, 


5 The notation 7,1 is avoided because the Np-th gap is followed by the first and not by an (Np+1)-th gap. 
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(a) Full-arc bearing with Np = 0 pads. (b) Multi-pad bearing with Np = 3 pads. 


= 


Figure 3.4: Sketches of exemplary domains on which the relevant fluid PDEs can be solved by referring 
to boundary conditions for inflow (solid lines), outflow (dotted lines), and symmetry (dot-dashed lines). 


the trailing edge of each domain coincides with the leading edge of the next domain, 
where the Np-th pad is followed again by the first pad. When splitting the circle T! 
in Eqs. (3.57), (3.58a), and (3.58b) accordingly, this leads to Np symmetry boundaries 


IOs = (Xn Xn+4x) x {0}, n=1,...,Np (3.60) 


for any multi-pad GFB and, analogously, to as many inflow and outflow boundaries 


Qni = ne) € (Xn, Xnt Ax) x 15} | - > o) U {Xn} x (6 3). (3.61a) 


dOn o = ne) € (Xn, Xn+AX) x 15} | r < o) U {xn+Ax} x (05). (3.61b) 


In Egs. (3.61a) and (3.61b), the additional inflow regions at g = Xn and the additional 
outflow regions at 9 = x, +Ax can be included without any case distinction since 
the velocity component ix in Eq. (3.33a) is always positive® for sufficiently large wo. 


3.4.3 Selection of Boundary Conditions 


On the inflow boundaries, the thermodynamic state of the fluid is supposed to be 
entirely known, which dictates the ambient density 0%, the ambient temperature do, 
and all other properties owing to thermal, caloric, and viscous equations of state. 
On the outflow boundaries, which are only relevant for parabolic terms, the ambient 


6 This supposes a positive rotational speed, but it is trivial to simply swap the edges in the contrary case. 
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pressure pæ alone acts back on the flow in the upstream direction, while density, 
temperature, and all other properties are extrapolated from the domain as required. 
Besides these Dirichlet boundary conditions, homogeneous Neumann boundary 
conditions are introduced on the symmetry edges and make disappear the gradients 
of all properties across the bearing centerline. Unsurprisingly, this entire approach 
shows many parallels to the usual procedure for dealing with general compressible 
flow problems (see, e.g., [FPS20]), which is transferred here to fluid film lubrication. 


While the specification of constant Dirichlet boundary conditions is common practice 
at the axial bearing edges (see, e.g., [Sze1o]), the validity of such an assumption may 
well be questioned considering the foil fixation gaps. This uncertainty is due to the 
possible mixing of fluid from the environment with the partially recirculating fluid 
in the lubrication gap, which would require an empirical mixing parameter that is 
difficult to determine and strongly related to the exact fixation gap geometry [Gar12]. 
However, knowing that a good thermal design should ensure the most complete fluid 
exchange possible according to San Andrés et al. [SK10], such ambient conditions at 
the fixation gaps are easily justifiable when thinking of idealized bearings anyway. 


Combining all of the previous assumptions and referring to domains and boundaries 
that are defined in Eqs. (3.56), (3.57), (3.58a), (3.58b), (3.59), (3.60), (3.61a), and (3.61b), 
the fluid BVPs to be solved can ultimately be stated for any number of pads as 


Eq. (3.53), in Oy 
Eqs. of state 


P = Po, 0 if Np =0, 
0 = eo, Os ( PETAS G ...,Np if Np > 1. 6.62) 
Egs. of state 


p = Po on 0On 6 
9/92 =0 on 0Qn,s 


Here, it should of course not be forgotten that the radial boundary conditions from 
Eqs. (3.23a), (3.23b), (3.29), and (3.47) are already incorporated into the respective 
PDEs as an outcome of them being integrated from one bounding wall to the other. 
From Eq. (3.62), the pressure distribution along the entire circumference as required 
for the calculation of bearing forces in Eq. (2.25) is obtained by sequencing all the 
domain-specific solutions together with the values of Dirichlet boundary conditions. 
To finally complete the problem formulation, suitable initial conditions are specified 
for the density and temperature fields, so that all other properties can be deduced 
via the respective equations of state. In subsequent investigations, this implies that 
the ambient state from 0(Q;,; in Eq. (3.62) initially prevails throughout the domains, 
which should of course not be mistaken for the description of an actual GFB start-up. 
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Dry Friction Characteristics 


As the most obvious feature, the compliant metal foil structure gives such flexibility 
to the outer fluid film boundary that GFBs can compensate manufacturing tolerances 
and may adapt to different operating conditions much better than rigid gas bearings. 
Under fluctuating pressure load, the foil strips are subject to minute relative motion 
accompanied by dry friction, thus providing an adaptive damping mechanism that 
reduces or even prevents undesirable vibrations of the entire bearing-rotor system. 
To reproduce the relevant nonlinear effects while staying computationally efficient, 
this chapter first reduces the number of parameters to a minimum before proposing a 
lumped-element bump foil model of reasonable complexity. Based on this approach, 
it is feasible to describe friction in various ways reaching from a regularization of 
Coulomb’s law of friction to dynamic bristle models capturing stick-slip transitions. 


4.1 Lumped-Element Modeling Approach 


Besides the lubricating fluid, the foil structure represents the main difficulty but also 
the main interest in modeling GFBs, which is mostly due to its particular design 
with a multitude of resulting effects. Since there is consequently no lack of modeling 
attempts in the literature, many of which seem overly complex and rather inefficient, 
it is all the more important to first clarify the relevant requirements. In this thesis, 
one of the objectives consists in a systematic understanding of how friction affects 
and possibly improves the overall system dynamics. Therefore, the primary purpose 
of the established minimal model is to qualitatively reproduce this strongly nonlinear 
influence while introducing the smallest possible number of empirical parameters. 
In doing so, given that friction is directly related to the occurrence of relative motion, 
it becomes of course essential to take into account not only local displacements of the 
compliant structure but also the far-reaching interaction between deforming bumps. 


The schematic sketch in Fig. 4.1 shows the considered lumped-element model, which 
approximates the basic geometry and thus also the kinematics of a bump foil strip by 
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Figure 4.1: Sketch of the lumped-element foil structure model for an exemplary first pad with five bumps, 
in which the bump foil (in light green) is reduced to an arrangement of rods, springs, and mass particles, 
whereas the top foil (in dark green) is imitated by simply interpolating linearly between adjacent bumps. 


representing each bump as a triangular arrangement of two rigid rods and a spring. 
As already explained in Section 2.2 in a more general context, this structure model 
is not only imagined “unwrapped” in conformity with the fluid model, but it is 
also a plane model that is sufficient for an axially almost constant deformation field. 
Referring to dimensions that are frequently found in the literature, the undeformed 
bump length* 2bg < RAY and the undeformed bump height hg define the geometry, 
which yields the rod length \/(b3 + h3) by consideration of an undeformed bump?. 
Idealized as roller supports in the sketch, the leading edge and the connecting rods 
between the triangles can move horizontally but cannot detach from the surface, 
whereas the trailing edge rests on a pinned support in imitation of the foil fixation. 
In this greatest possible simplicity, which nevertheless allows to include friction and 
bump-bump interaction, the described model is similar to that of Feng et al. [FK10]. 


The vertical forces acting on the bump apices are supposed to represent the resulting 
pressure forces that, in reality, are transmitted via areas in contact with the top foil. 
Depending on the local differences between the pressure in the lubrication gap and 
the assumed ambient pressure between and below the foils, the top foil adapts to the 
deformed shape of the bumps and thereby sags into the clear space between them. 
In a preliminary work by the author [Lei15], the aforementioned foil sagging effect 


* It may surprise in view of Eq. (2.4) that the undeformed bump length is not simply equated with RAY. 
However, as can be seen in the course of the nondimensionalization in Section 5.1, the new quantity 2b, 
is actually introduced on the radial length scale C rather than on the circumferential length scale R. 
Hence, both lengths are independent but should fulfill the given compatibility condition for consistency. 

? Note that the fourth bump from the left in Fig. 4.1 is actually shown undeformed because of x14 = x15. 

3 It goes without saying that the model discussed here is not suitable for describing full bearings without 
any fixation gap, which makes the limiting zero-pad case Np = 0 only applicable with a rigid structure. 
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is thoroughly investigated by modeling the top foil as a segmented beam structure. 
The found results show that the wavy shape of the deformed top foil is reflected in a 
waviness of the fluid pressure field even though this has little or no influence on the 
qualitative rotordynamic behavior. When discretizing the lubrication gap, however, 
the local character of this effect requires an exaggeratedly fine computational grid, 
which seems hardly justifiable compared to the lumped-element bump foil model. 
In the present thesis, the top foil segments are imagined as being extensible rods that 
straightly connect bump apices and thereby represent the outer fluid film boundary. 
For the sake of completeness, Appendix A shows the effect of top foil deformation on 
the results from Section 6.3.1 and thus builds a bridge between this work and [Lei15]. 


From an engineering perspective, lumped-element modeling of continuous structures 
is only permissible if no natural frequencies near the expected excitation frequencies 
are lost (see, e.g., [DF20]). According to numerical and experimental results reported 
by Bin Hassan et al. [BB17], the first natural frequency of a typical bump foil is found 
in the region of 2 kHz, that of an isolated single bump is even about ten times higher. 
Given that the frequency range of interest correlates with the rotor speed and thus is 
comparatively low, the foil mass seems to be completely negligible at first glance. 
However, practical experience suggests that uncertainties involved in the treatment 
of DAEs are far more delicate than challenges arising from solving systems of stiff 
but ordinary differential equations [OLS17]. Hence, the roughly approximated mass 
of individual bumps may be concentrated in discrete mass particles with 


mp = pprtbgogL. (4-1) 


The exact value of Eq. (4.1) being insignificant, the foil density og is multiplied by the 
bump-like volume of a half-cylinder shell of radius bg, thickness og, and height L, 
where 03 is interpreted as the foil thickness and L denotes the bearing length. 


The search for an appropriate definition of the spring stiffness is also limited to a 
rough estimation of an order of magnitude, thus providing a range for parameter 
variations without restricting the analysis to a particular foil design and geometry. 
Based on the pioneering approach by Walowit et al. [WA75], which is still widely 
used despite its surprising simplicity and the many suggestions for improvement, 
each bump is modeled as an elastic circular arc under axial plane-strain conditions. 
Depending on Young’s modulus Eg and on Poisson’s ratio vg, the stiffness in relation 
to a force applied at the apex is then given in first approximation by the expression 


EgL OB ? 
k = — . . 
B ae) (4.2) 


In fact, the simplification behind Eq. (4.2) is essentially based on the finding that the 
curvature of bumps can be neglected without introducing substantial inaccuracies. 
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Table 4.1: Properties of a typical foil structure. 


Description Symbol Value Unit 
Mass density PB 8.28 g/ cm? 
Young’s modulus Eg 214 GPa 
Poisson’s ratio VB 0.29 — 

Foil thickness OB 0.10 mm 
Bump half-length bg 178 mm 
Bump height hg 0.51 mm 
Bearing length L 38.10 mm 


The inclusion of nonlinear bump stiffness formulations as proposed by Baum [Bau18] 
is straightforward and may improve quantitative results for practical applications. 


Concerning realistic dimensions of the bearing and of the foil structure, the report 
by Ruscitto et al. [RMG78] is among the first to give an overview of such properties. 
Ever since, this data is of regular use in the literature* and is understood as reference 
for the comparison of results, so that no deviation from this practice is intended here. 
The report mentions the nickel-chromium INCONEL alloy X-750 as the foil material, 
which finally yields Table 4.1 when referring to the manufacturer data sheet [Speo4]. 
Using these parameter values, the bump mass mp equals a few tenths of a gram and 
the bump stiffness kg is around 1 N/m, which is both in agreement with [BB17]. 
Lastly, the rod length is only slightly larger than bg for the small bump height hg, 
meaning that this geometry parameter is expected in the region of a few millimeters. 


4.2 Foil Strip Kinematics and Interpolation 


Seeing that the lumped-element model of the n-th pad? has Ng degrees of freedom, 
the configuration of a given foil strip can be described by Ng generalized coordinates. 
Here, the circumferential displacements® x,,1 of the mass particles shown in Fig. 4.1 
seem the most natural choice for a minimum number of independent coordinates. 
Compared to other selections, each generalized coordinate xy, with the resulting 
generalized velocity %nm later allows for an immediate characterization of friction 
between the bump foil region that precedes the m-th bump and the bearing sleeve. 


4 The original dimensions being reported in inches, values in SI units differ slightly between publications. 
In the present thesis, the converted properties listed in Table 4.1 are also rounded as deemed reasonable. 


5 Since all pads are modeled identically, it is sufficient to consider an arbitrary n-th pad in isolation. 
Consequently, the index n is always mentioned hereafter but it remains unchanged until further notice. 


6 Note that Xn,m is introduced in the negative e,-direction, which later reverses the signs of friction forces. 
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Moreover, even though not explicitly needed for the described modeling approach, 
the circumferential displacements and velocities of bump apices may be obtained by 
just averaging the two concerned generalized coordinates or velocities, respectively. 


As indicated by the sketch in Fig. 4.1, the radial displacements yn,m of the bump 
apices are expressed relative to their respective positions when all springs are relaxed. 
Using the Pythagorean theorem to determine the deformed height of the m-th bump, 
the apex displacement depends on generalized coordinates by a kinematic relation 


Ynm = Jala H 3) (26 F Xn,m nmi) = hp: (4-3) 


To avoid case distinctions in Eq. (4.3) and other equations to come, the symbol x, n3+1 
is an alias for the vanishing displacement at the foil fixation and must be set to zero. 
Requiring a mapping at the differential level, it reveals useful to define the Jacobian 


OY nl sees Yan J nA =, J nl 0 
OXn,1 IXn,Np In2 —Jn2 
= : at : = In3 —In3 j (4.4) 
OVn,Ng ae OVn,Ng ) is 
OXn,1 OXn,Ng 0 Ja, Ng 


which corresponds to an upper bidiagonal matrix because Eq. (4.3) involves only the 
generalized coordinates xy, and xX, +41 of the two directly adjacent mass particles. 
Moreover, given that Eq. (4.3) is anti-symmetric in the two generalized coordinates, 
the matrix in Eq. (4.4) is entirely determined in terms of only Ng diagonal entries 


_ OYn,m u 1 2bg + Xn,m — Xn,m+1 


Dm * Jan) (2b5 + nm innit) 


Taking finally advantage of Eq. (4.4) when analyzing foil motion, the transformation 


(4.5) 


Inm 


Un Xn 
Yn2 Xn2 
Yn3 | = In Xn3 (4.6) 
Yn,Ng Xn,Np 


yields all the radial velocities 1, of the bump apices, which depend not only on 
the generalized velocities, but through Eq. (4.5) also on the coordinates themselves. 


Geometrically speaking, the top foil corresponds to straight lines that permanently 
connect the bump apices, so that position and velocity of any point can be obtained 
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through linear interpolation based on yn,m and nm of the two concerned bumps. 
At fixation gaps between neighboring pad edges, it seems fair to extrapolate values 
of the closer bump apex, thus creating a regularizable discontinuity in the middle. 
With Eqs. (4.3) and (4.6), this makes the deformation and deformation velocity fields 


a d 
(pt) = qe. Xnm.), (Pt) em) (47) 


functions of the generalized coordinates and velocities, the piecewise definitions of 
which are not explicitly stated above because they are of purely technical nature’. 


4.3 Derivation of Equations of Motion 


The established foil structure model corresponds to Np systems of Ng mass particles, 
for which Np systems of Ng coupled equations of motion need to be formulated. 
With the numerous constraints and the repetitive geometry, it seems advantageous 
to proceed for the n-th pad via the Euler-Lagrange equations (see, e.g., [See18]) 


d f In OL, 
= =; : 
dt & | ( | Qnm i posg Ny 48) 


which involve the Lagrangian £n and also a number Ng of generalized forces Qn m. 
The kinetic and the potential energies for the n-th pad being particularly simple, 


it comes 
Ng Ng 2 
mpg 2 kg 
Ln = ar L Anm — > (re = hel) , (4.9) 
m=1 m=1 


an expression in which gravity does not appear for the “unwrapped” structure model. 
Moreover, formulating the virtual work of the non-conservative fluid forces Fy im and 
friction forces fn,m in generalized coordinates, a comparison of coefficients yields the 
generalized forces 


On Fa fna 
On2 F2 fn2 
Qn | =J} | Fna | + | fn |. (4.10) 
Qn, Ng Fn, Ng Fn,Ng 


7 The interested reader is provided with detailed formulae of such a linear top foil interpolation in [Bar16], 
which represents a preliminary work inspired and co-supervised by the author of the present thesis. 
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Most notably, Eq. (4.10) takes advantage of the transposed Jacobian J], from Eq. (4.4) 
to relate pressure forces that act on the bump apices to the generalized coordinates. 


Inserting Eqs. (4.9) and (4.10) into Eq. (4.8) and calculating the respective derivatives, 
the obtained outcome for the n-th pad is the system of coupled equations of motion 


Xn 1 -1 0 Xn,1 
Xn2 —1 2 —1 Xn,2 
Ën, + we —1 2 -1 Xn,3 
Ën, Ng 0 —1 2 Xn, Ng 
T 
Jna Jna 0 Faa fna 
i Jn2 —Jn2 Fn2 i fn2 
Ss In3 Ina Fa3 — — Ín j (4.11) 
MB p mB 
0 Jn, Ng Fa, Ng Ín, Ng 


in which the inherent dynamics are characterized by the natural angular frequency 


wB = er (4.12) 


While the linear left-hand side of Eq. (4.11) represents a harmonic oscillator chain, 
the right-hand side corresponds to nonlinear fluid forces and frictional damping. 


To assess the resulting fluid forces that are exerted on the respective bump apices, 
the locally load-distributing function of a top foil is roughly imitated by integrating 
elementary pressure forces over the associated area of negligible curvature. With the 
ambient pressure pœ underneath the foil structure counteracting the pressure field p 
in the lubrication gap, this becomes 


E rtm + Ap 
Fe 2f / (p — po) Rdg dz. (4.13) 


Since any circumferential displacements of the foil structure are judged negligible 
compared to the length scale Zo of overall bearing dimensions, the integration limits 
in Eq. (4.13) rely on Eqs. (2.4) and (2.5) that refer to the undeformed foil geometry. 
Eventually, considering the fluid forces from Eq. (4.13) for the solution of Eq. (4.11) 
and inserting the kinematic relationships from Eq. (4.7) into Eqs. (2.13) and (2.15), 
the occurring fluid—structure interaction is entirely captured in both directions. 
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4.4 Discussion of Dry Friction Formulations 


As a fundamental design feature of GFBs, which are prone to undesirable vibrations 
due to the low viscosity of gases, the compliant and movable foil structure has the 
beneficial property of introducing dry friction as additional damping to the system. 
For the description of friction-induced effects, which is of course not only a matter 
in GFBs but also in many other mechanical systems, the available literature offers a 
very large selection of models in varying levels of complexity for different purposes. 
A general overview of relevant approaches is given in the thesis of Kapelke [Kap19], 
while a comprehensive review of many friction models can be found in [Mar+16]. 


Relative motion and therefore dry friction occur in bump-type GFBs principally in 
two different types of contact regions, which are located between bearing sleeve and 
bump foil on the one hand and between bump foil and top foil on the other hand. 
In connection with the present work, a preliminary finite-element analysis [Alv18] 
distinguishes exactly between the respective influences of friction at this or that spot, 
which leads to the simple conclusion that all observed effects are practically the same. 
Consequently, and also in the very spirit of a lumped-element modeling approach, 
the only concentrated friction forces fm in Fig. 4.1 act directly on the mass particles. 
Setting the normal forces equal to the nearby pressure forces F;,1 from Eq. (4.13), 
a combined expression for Coulomb’s law of dynamic friction (with coefficient jg) 
and for a corresponding regime of static friction can be given in the set-valued form 


{—upEnm } if Xn,m < 0, 
fnm € en: +p Fam] if tnm = 0, (4.14) 
{+upEnm } if Xn,m > 0. 


In Eq. (4.14) and also in other upcoming friction formulations, the normal force Frm 
is assumed to be positive, so that the concerned friction force is set to zero otherwise. 
As a matter of fact, this can be considered a basic representation of foil detachment, 
knowing that more sophisticated contact mechanics approaches [AB19] are difficult 
to be integrated into a model that ultimately results in a single dynamical system. 


To further complicate the problem, the insertion of Eq. (4.14) into Eq. (4.11) leads to 
a system of differential inclusions, which can roughly be seen as multi-valued ODEs. 
For such Filippov-type equations, especially the study of bifurcations is an own field 
of research that is thoroughly discussed, for instance, in the thesis of Leine [Leioo]. 
Since the GFB-rotor model in the present thesis reveals far too complicated to be 
treated using such a method, it is discarded in favor of some regularized formulation. 
In this context, it is common practice to neglect static friction and to express Eq. (4.14) 
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by a sign function, which can then be approximated by a hyperbolic tangent function. 
Introducing a small regularization velocity ug < Xn,m, this transforms Eq. (4.14) into 


x 
Fam = UBEnm tanh (=) (4.15) 


and can thus be included in Eq. (4.11) without affecting the solution requirements. 
On the downside, it must of course not be forgotten that this convenient formulation 
has some deficits in view of possible sticking conditions, however, it still represents 
an important model improvement when compared to equivalent viscous damping. 


For the sake of completeness, it must be mentioned that more advanced formulations 
from the class of bristle models can of course also be considered, including notably 
the Dahl model [Dah76], LuGre model [Can+95], and elasto-plastic model [Dup+o2]. 
The idea behind such approaches basically consists in the introduction of a tangential 
contact compliance to account for pre-sliding deformations of the surface asperities, 
which thereby provides some kind of physically motivated regularization [Kap1g]. 
Applied to the foil structure model, each bump needs an additional coordinate® zn m 
to describe the linearly elastic deflection (with stiffness kz) of the associated bristle. 
Assuming for brevity a most simple evolution rule for this internal bristle variable 
following the Dahl model [Dah76], the friction equations coupled to Eq. (4.11) write 


Tnm 
uBEn,m 


Jam = KzZnms Znm = Xnm — (tnm (4.16) 
Especially when the right-hand side of the ODE in Eq. (4.16) is modified by adding 
a suitable switch function as in the elasto-plastic model by Dupont et al. [Dup+02], 
such a friction model reveals capable of actually reproducing stick-slip transitions. 


Nevertheless, no model that truly accounts for sticking can circumvent the inherent 
and physically originated difficulties that inevitably arise when facing dry friction. 
Even though Eq. (4.16) does not lead to a differential inclusion as does Eq. (4.14), 
another directly related issue is observed upon consideration of the evolution rule. 
More precisely, the previously mentioned switch function in the elasto-plastic model 
must be defined such that Eq. (4.16) becomes żn,m = Xn,m during phases of sticking. 
In this case, the concerned bristle effectively loses its additional degree of freedom, 
so that suitable solution strategies have to deal with varying problem dimensions. 
Since it is not intended to perform stability and bifurcation analyses for such systems 
within the scope of this thesis, most of the presented findings are based on Eq. (4.15), 
while some time integration results with a bristle model are also published in [LSB18]. 


8 The symbol Zn,m is chosen in reference to the literature and has nothing to do with the axial direction. 
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5 Computational Analysis 
of GFB-Rotor Systems 


The comprehensive fluid-structure-rotor modeling approach that is described in 
previous chapters leads to a multi-dimensional system of coupled PDEs and ODEs, 
so that an analytical treatment does not seem to be a very promising undertaking. 
Consequently, this chapter outlines the way toward a numerical solution strategy, 
according to which the equations are nondimensionalized and spatially discretized 
to give a mathematical dynamical system that can be investigated computationally. 


5.1 Nondimensional Problem Formulation 


Owing to the coupled nature of the problem, it appears indispensable to establish a 
nondimensional formulation that remains consistent throughout all three submodels. 
First of all, the introduction of a characteristic time scale to yet to be determined 
defines the nondimensional time 


T= iy’ (5.1) 


with the corresponding nondimensional time derivative denoted by (...)’ for short. 
Coming to the independent Cartesian coordinates that describe the lubrication gap, 
a suitable nondimensional circumferential coordinate ¢ is already given in Eq. (2.6) 
and is completed, with the bearing length L, by a nondimensional axial coordinate 


Z 
Z = —., . 
L (5.2) 
In the direction of the film thickness, the radial clearance C is the natural reference 
to define nondimensional forms of all concerned kinematic quantities in Chapter 2, 
most notably foil structure deformations and rotor displacements, according to 
h q 1 e 
H= =, NE X=, Y=, NS =T. C 
C = C C TE y (5.3) 
With Eq. (5.3), any dependencies on the dimensional coordinates are also transformed 
to their nondimensional counterparts, which implies y(t) = T (Tto) for instance. 
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For the nondimensionalization of fluid BVPs, the relevant thermodynamic properties 
in Chapter 3 are related to the respective values under normal conditions by setting" 


p p v u w 
P=—, D=—, T=, V=—, W= —. (5.4) 
5 2 5 m D 5-4 


Then, it seems a natural choice to express the nondimensional derivative properties 


°9° 
Ay = day, Cv = z 


Cy (5-5) 


using the same reference values as in Eq. (5.4) and combining them appropriately. 
Even though not stated explicitly for the sake of brevity, such nondimensionalizations 
apply analogously to fluid boundary conditions such as To = Bo/0° and others. 
Taking advantage of Eqs. (5.1) to (5.5), the PDE from Eq. (3.36) results immediately 
in the nondimensional Reynolds equation for compressible fluids 


9 Ag à _ ð (DH?OP\  , 0 ðf DH oP 
a (DH) +5 ap (PR) = 2V 3 a waz} 59 


just like the PDE from Eq. (3.53) yields the nondimensional temperature equation 


DC,H 


dr 209 2W*i\dpdp | dZdZ 
ƏH , Ao0H 9 (HaP\ 20 { H°0P 
OT 2 dp dp|\2V dp 09Z \ 2V ƏZ 


oP : oP ; 
GEO 


Having a closer look at Eqs. (5.6) and (5.7), it becomes clear that their particular form 
fixes not only the time scale 
6u’ (R 3 
to = =], 8 
= F(z) 68) 
but also results in an appearance of the three nondimensional similarity parameters 


R 6u° (R\2 ve Rv 
KO = 7 Ag = towo = = (à) Wo, © = le) Bo- (5-9) 


aT  AoûT (er 3) 


= -A,TP 


AV, HB 
6H 22V 


&(T- Ta) + 


* This supposes the vapor quality to be larger than 0 % under normal conditions, which holds for R-245fa. 
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In Eq. (5.9), the first parameter xo is a purely geometric number referring to the ratio 
between radius and length of a GFB, which distinguishes short from long bearings. 
The second parameter Ag resembles the Sommerfeld number (see, e.g., [Sze1o]) and 
is subsequently denoted as bearing number or nondimensional rotational speed. 
As an important matter of fact, the bearing number is the only place in the overall 
nondimensional problem in which the angular frequency wo of the rotor is involved, 
while in particular the adopted time scale in Eq. (5.8) is independent of the speed. 
For run-up and coast-down simulations or when considering Ag the bifurcation 
parameter in solution continuation, this ensures to not misleadingly bias the results. 
Finally, the third parameter ©, in Eq. (5.9) is called thermal management number 
and describes the general ability of a GFB to transfer heat to the environment. 


Concerning the foil structure model, the nondimensional form of displacements 
should be compatible with the nondimensional deformation field Q from Eq. (5.3), 
so the clearance C is used in the circumferential and in the radial direction to give 
= #8 


Hg = TC (5.10) 


Xn,m Yn,m Zn,m 
Xn,m = , Yım = , Zum = , Bg = 


C C 


Any loading of the foil model being related to surface forces due to pressure, 
nondimensional forms of Eqs. (4.13) to (4.16) are given by the respective expressions 


Fam m 
Hya = Jom 
nm RL mm RL 


(5-11) 


Altogether, this transforms the system of equations of motion from Eq. (4.11) into? 


Xna ] | Xn,1 | | Ina ] | Tn, ] 
n,2 Xn,2 T2 Tn,2 


XF 3 + (en |; 5 | Xn3 = — Me l. 7 | [3 = Me TTn,3 A (5.12) 
5 R B . B 3 
aN Xn, Ng Iln Ng Tn, Ng 


so that the foil behavior is characterized in terms of two nondimensional parameters 


6u? R 2 pe C 5 
ete Zu ; 
B = towg D (à) Wp, Mpg 36G JL (à) mp, (5.13) 


referred to as nondimensional bump frequency number and bump mass number. 
Besides Eq. (5.13), the nondimensional friction coefficient ug in Eqs. (4.14) to (4.16) 


? Since the two large matrices in Eq. (4.11) are formally identical to their nondimensional equivalents, 
they are not repeated here. In the same way, the large matrix in Eq. (2.23) is later omitted for the rotor. 
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represents a native nondimensional number that is unaffected by any reformulation. 
When considering a bristle friction model, however, the ODEs from Eq. (4.16) become 


TCn,m 


uBlIn,m | 


1 2 2 
Tnm = KzZnm Zim = Xn,m — bean 


(5.14) 


in which the definitions of nondimensional displacements and forces in Eqs. (5.10) 
and (5.11) make the bristle stiffness number Kz = C/(p°RL)kz another parameter. 


Regarding the rotor model, the same nondimensionalization as in Eq. (5.3) leads to 


_ SR _ MR _ ER 
XR = C” Yr Cc’ ER= E (5-15) 


and the same consideration as in Eq. (5.11) gives the nondimensional bearing forces 


Hy — 6 = 
x PRL’ YS 


Fy 6 
p°RL x (5.1 ) 


Applying this to Eq. (2.23) results in a system of nondimensional equations of motion 


x” 0 X 
y"|  2Dr lol -O8 | | Y 
X#|  1-Er [X| 1-2r l] | Xp 
xs Yk YR 
IIx 0 0 
2 IT, 2 0 1 
2e — ERA — : 
eme | 0 | 849 |sin(agz)| lol 647 
0 cos(AgT) 1 
with four newly appearing nondimensional numbers 
© 2 (0) 2 
OR = towR = (à) WR, Dr = toûr = Sr (à) ÖR, (5-18a) 
o 5 o\2 5 
p C 36(u°) [R 
= = = = x . 8b 
Me = zur) Pr TA TE sine 


While OR and MR are constructed and interpreted in the same way as Qg and Mg 
from Eq. (5.13), the characterization of the rotor requires two additional parameters. 
To this purpose, the viscous damping is reflected in DR according to Eq. (5.18a), 
while G in Eq. (5.18b) represents the gravitational loading of a horizontal rotor. 
Obviously, the discussed nondimensionalization applies analogously to Eq. (2.24) 
for the rigid rotor model, which is then independent of the frequency number Qpr. 
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5.2 Finite-Difference Discretization Schemes 


While the foil structure and the rotor are represented by discrete models that lead 
directly to systems of exclusively time-dependent ODEs, the continuous fluid model 
produces PDEs that depend on both time and space. Lacking analytical solutions, 
the common remedy in computational fluid dynamics (CFD) consists in using spatial 
discretization schemes to transform such PDEs into systems of coupled ODEs that 
deliver approximate solutions. Among the different classes of available approaches, 
e.g., finite-element (FEM), finite-volume (FVM), or finite-difference (FDM) methods, 
the selected FDM provides the greatest flexibility when facing non-standard PDEs. 


Since the nondimensional domains that describe the lubrication gap are rectangular, 

the introduction of a corresponding structured computational grid is not complicated 

when specifying the number of grid points in the respective directions by Ng and Nz. 

However, referring to the overall fluid region to define the equidistant grid spacings 
27 1 


A= NT AZ = x: (5-19) 


some additional care is indicated to ensure compatibility with the number of pads Np, 
with the number of bumps Np, and with the axial symmetry boundary conditions. 
More specifically, requiring grid lines that coincide exactly with each fixation gap, 
with each foil-sleeve contact line, and with the bearing’s circumferential centerline, 
a closer look at the exemplary grids in Fig. 5.1 yields two discretization conditions? 
Nz-1 


Ap _Nọ-1 3 _ 
Ag = WN; EN, WS EN. (5.20) 


Finally, having established the computational grid according to Eqs. (5.19) and (5.20), 
the grid origin is anchored to the nondimensional physical domain in such a way 
that the integer grid coordinates (i, j) correspond to discretized physical coordinates 


gi = (i-1Aptx, i=1,..., Nọ, (5.21a) 


Zj = (j-1)AZ 


1, 
5’ j=l,...,Nz. (5.21b) 
Inserting Eq. (5.21a) into the nondimensional forms of Eqs. (2.13), (2.15), and (4.7), 
these field quantities transform into a discretized fluid film thickness H; and into a 
discretized foil structure deformation Q; with discretized time derivatives H! and Q’. 


3 Graphically, this implies that Ay is an integer multiple of Ag, while 1/2 is an integer multiple of AZ. 
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Boundaries FDM stencil 


g=0 Symmetry 


(a) Full-arc bearing with Np = 0 pads (described (b) Multi-pad bearing with Np = 3 pads (described 
by the total number of Ngpm = 27 internal nodes). by the total number of Nepm = 18 internal nodes). 


Figure 5.1: Sketches of exemplary computational grids for the FDM with Ng x Nz = 10 x 7 grid points, 
which are either internal nodes (circles), boundary nodes (squares), or mirror nodes (dashed circles). 


Considering once more the sketched grids in Fig. 5.1, it becomes obvious that only 


New = T (Ng ae Np) (Nz _ 1) (5.22) 


grid points are internal nodes that lie in domains on which PDEs need to be solved. 
At each internal node, a discretized density Dj; and a discretized temperature T; ; 
are introduced as variables, on which all other discretized thermodynamic properties 
such as the pressure P; ; or the viscosity V;; depend via the given equations of state. 
Around the internal nodes, the remaining grid points are either boundary nodes, 
which store constant values given by Dirichlet boundary conditions, or mirror nodes, 
in which values from nodes of the other symmetric half of the domain are reflected. 


The basic idea behind any FDM consists in the replacement of partial derivatives by 
corresponding difference quotients, which is close to the definition of a derivative 
and can be justified by reformulating a truncated Taylor series (see, e.g., [FPS20]). 
Evaluated on the introduced computational grid, this makes any spatial derivative 
at the position of a given internal node (i, j) an algebraic function of discrete values 
obtained from the node itself and from four neighboring nodes (i+ 1,/) and (i,j +1). 
Regarding the PDEs in Eqs. (5.6) and (5.7), the diverse nature of occurring derivatives 
as discussed in Section 3.4.1 leads to a distinction between terms of parabolic and 
terms of hyperbolic character, which affects the choice of suitable approximations. 
Following an approach that is typically used in CFD when dealing with any kind of 
diffusion equation, the second-order terms with inner and outer derivatives on the 
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right-hand side of Eq. (5.6) can be expressed by central differences (see, e.g., [FPS20]) 


Din jH y Di- jH} y 
o (DH3aP\| _ Win (Pau, P) ARTE 2 (Pij PA) Bi 
ae | 2V dg }\ 7 Ag? FE" 
Lj 
D; 1, H? Di; H? 
tat; ij—V2 ti #2 > 
a (DH’OP\|_, In (Pin Bi) PAPE (P; E a 
az | 2v əz || ~ AZ? rn 
Lj 


Such schemes, in which the fractional indices i + 1⁄2 or j + Y2 denote mean values of 
adjacent nodes, are also used for the terms in Eq. (5.7) without density as a factor. 
For the first-order derivatives in Eqs. (5.6) and (5.7), from which the derivatives of 
temperature are picked out as an illustrative example, numerical stability becomes a 
relevant issue that can be addressed by using upwind differences (see, e.g., [FPS20]) 


Tija Tij à 
oT 5 Tij zei Ti, oT à A. if Pi j+1 > P;j-1 (5.24) 
dp ij Ay ðZ ij Hz otherwise. 


In order to determine the direction of flow in Eq. (5.24), the same consideration as for 
the distinction between inflow and outflow edges in Eqs. (3.58) and (3.61) is applied. 
Accordingly, the circumferential velocity is assumed to be always positive, whereas 
the axial velocity must be judged by means of the discretized axial pressure gradient. 
Of course, the first-order derivatives in Eqs. (5.6) and (5.7) that involve P, D, or H 
are discretized in analogy to Eq. (5.24) even though not stated for the sake of brevity. 


When the finite-difference stencil moves around to each of the Nppm internal nodes, 
the FDM-discretized Reynolds equation and temperature equation produce together 
a total number of 2Nrpm first-order ODEs, which are reordered to the abstract forms 

Di, = FDMp Di Di+1,j, Di jt, Ti 


Tizi, Ti 1 Hi, Hiss, H Hu}, (5-25a) 


Lj jr j 


14T: 


Tij = FDMr{ D;;, Dji41,j Dij jr Tizi, Ti j1 Hi Hizi, Hi, Hi} (5.25b) 


Realizing that peripheral nodes of the stencil may fall on boundary or mirror nodes, 
it becomes clear that these equations represent not only PDEs but even entire BVPs. 
Most remarkably, it is finally the very nature of the selected schemes that enables 
the central-difference parabolic terms to consider boundary conditions at all edges, 
while the upwind-difference hyperbolic terms do not even “see” outflow boundaries. 
Consequently, any solutions of the discretized ODEs in Eq. (5.25) are also considered 
approximate solutions of the original BVPs from Eq. (3.62) in nondimensional form. 
In order to assess the quality of the approximation, grid independence studies as 
shown exemplarily in Appendix B are conducted prior to performing actual analyses. 
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5.3 Nonlinear Dynamical System Representation 


Although the establishment of a dynamical system representation is a decisive step, 
the feasibility of which is an important merit of the present thesis, this operation does 
not cause major difficulties in view of the preparatory work already accomplished. 
It only remains to identify suitable variables that uniquely describe the system state 
when forming a smallest possible fluid-structure-rotor state vector (or phase vector) 


es 
s=| $ $ sk |enck, (5.26) 


where the associated phase space (or state space) Q is an N-dimensional manifold*. 

Together with some evolution rule that characterizes the phase velocity according to 
dS N 

T = R{S, A}, R:QxR> R“, (5.27) 

this creates an autonomous dynamical system with the right-hand side vector field R 

depending solely on the current state and parameterized by the bearing number A. 


The fluid description? goes back to Section 3.1.1, where density and temperature 
are selected as independent primary variables on which all other properties depend. 
After nondimensionalization and discretization, this results in the fluid state vector 


T Nọ—1-Np)(Nz—1 
E ee eee € Rx = Rl NZ), (5.28) 


for which the evolution through time is known from Eq. (5.25) by first-order ODEs. 
In the case of an isothermal calculation that disregards the temperature equation, 
the discretized temperature is omitted and the fluid state vector is only half as large. 


With regard to the equations of motion that describe the foil structure in Eq. (5.12), 
each of the second-order ODEs generates one position and one velocity state variable, 
to which the respective bristle variable is added to give the foil structure state vector 


s 
ss= | es Re K Za tees, eR (5.29) 


2 r 


Here, the non-trivial evolutions through time of the second and third state variable 
for each bump can be obtained immediately from Eqs. (5.12) and (5.14), respectively. 
When the regularized Coulomb friction model is used instead of a bristle model, 
there is of course no need to include bristle variables in the foil structure state vector. 


4 The restriction of S to a subset of IRN relates to the requirement of non-negative Di; and T; in Eq. (5.28). 
5 Due to the symmetry of the problem, it is sufficient to consider fluid and structure of only one bearing. 
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Coming to the equations of motion of the rotor in Eq. (5.17), it is once more indicated 
to introduce position and velocity state variables for each of the second-order ODEs. 
Furthermore, it is usually purposeful and therefore required by Eq. (5.27) that all the 
ODEs are autonomized, so the resulting dynamical system is no longer explicitly 
time-dependent via the unbalance. Following a suggestion by Doedel et al. [Doe+o2], 
this can easily be achieved by the introduction of two artificial nonlinear oscillators 


Osin = +A00cos — Osin (in PO 1), (5.30a) 

Olos = —AgOsin — Ocos (Odin + Ožas — 1) (5.30b) 

and by replacing harmonic forcing in Eq. (5.17) by the asymptotically stable solutions 
Osin = sin(AgT), Ocos = cos(AoT). (5.31) 


Altogether, this yields the rotor state vector 
1 1 ! 1 d: 10 
Sr = | XoXo VV Xe KL Ye YO Où 0 | ER, (5.32) 


of which the evolution through time is governed by Eqs. (5.17) and (5.30), respectively. 
If no unbalance is considered, obviously no autonomizing oscillators are required, 
just as there is no need for disk positions or velocities when dealing with rigid rotors. 


Even though the most general form of Eq. (5.27) suggests that each component of R 
may depend explicitly on all of the N state variables, this is obviously not the case 
regarding the underlying equations. On the contrary, relevant coupling mechanisms 
can be broken down to the relatively simple interaction chart depicted in Fig. 5.2, 
which clarifies that there is no direct relation between the rotor and the foil structure. 
Most remarkably, the identification of physical quantities to be exchanged across 
domain interfaces allows for a modular design of the implemented research code, 
so that three submodels can be merged arbitrarily into one coupled overall model. 
Besides the most sophisticated modeling approaches on which this work is focused, 
the established collection of fluid-structure-rotor submodels includes, for instance, 
(visco-)elastic foundation models [LBS17b] and clamped rotor shaft models [LBS18]. 
Moreover, it is strongly emphasized that the described modular concept refers only 
to the assembly of governing equations, which are then solved simultaneously using 
a rigorous monolithic approach not to be confused with partitioned co-simulation. 
In this context, a recent contribution by Bou-Said et al. [Bou+20] gives a comparison 
of the two basic methods, which outlines the advantages of monolithic approaches 
in terms of precision and numerical stability, but also mentions possible limitations. 
Especially for less abstracted GFB geometries or even more complex fluid behavior, 
there would practically be no alternative but to choose a partitioned solution strategy. 
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Shaft Journal Kinematics Top Foil Kinematics 


Q(p,T), Q'(p,T) 


Fluid Submodel 
(Ideal Gas/ 


R-245fa/...) 


Bearing Force Components Distributed Fluid Forces 


Figure 5.2: Interaction chart illustrating the modular concept of the multiphysics model and underlying 
coupling mechanisms between the fluid submodel, the foil structure submodel, and the rotor submodel. 


The present thesis is not intended to provide an exhaustive documentation of the 
implemented research code GFB++, but here are some insights into fundamental 
technical aspects and underlying ideas. Basically, there are two project-specific main 
requirements that affect the choice of programming languages and environments, 
which are a sufficient computational speed to deal with the relatively complex model 
and an object-oriented programming paradigm to facilitate a modular conception. 
Therefore, the core code that evaluates the model equations is programmed in C/C++ 
and completed by Unix shell scripts for sequential program runs and post-processing. 
Essentially, when inheriting fluid-structure-rotor submodels from abstract classes 
that principally ensure the accessibility of all coupling quantities according to Fig. 5.2, 
solution algorithms can treat them as black boxes without knowing internal details. 
Finally, some relevant program parts are actually based on AUTO 2000 [Doe+02], 
which is a widely used package for solution continuation and bifurcation analysis, 
as well as on the C++ libraries Eigen [GJ+10], odeint [AM11], and CoolProp [Bel+14]. 
On a standard desktop computer, the expected computational times typically range 
from several minutes up to only a few hours for comprehensive parameter studies. 
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With the mathematical formalization of allmodel equations as one dynamical system, 
the associated theory provides a framework of analytical and numerical methods 
that are very extensively covered in the literature (see, e.g., [GH83, Kuzo4, MV11]). 
As a starting point for analyzing GFB-rotor systems, a nominal set of parameters 
is established before investigating the general stability and bifurcation behavior. 
Hereafter, particular features of refrigerant lubrication and also some possibilities 
of vibration reduction are thoroughly examined in view of the overall performance. 
In doing so, this closing chapter not only demonstrates the merits of the adopted 
modeling approach, solution strategy, and implementation, but also guides through 
a representative selection of the most important nonlinear effects and phenomena. 


6.1 Specification of Nominal Parameter Values 


The fluid-structure-rotor model and therefore the dynamical system in Eq. (5.27) 
depend on a large number of nondimensional parameters, for which physically 
meaningful values can be deduced from the corresponding dimensional parameters. 
As already mentioned at the end of Section 4.1 with regard to the properties of the 
foil structure, the present thesis is loosely based on GFB specifications published by 
Ruscitto et al. [RMG78], which are transferred here unchanged to three-pad bearings. 
For parameters with greater uncertainties, such as the heat transfer coefficient or the 
friction coefficient, moderate values that can reasonably be varied in both directions 
are chosen with some intuition. As further starting points for parameter variations, 
the ambient state corresponds initially to normal conditions and the standard rotor 
is as generic as those in preceding rotordynamic studies of Baum [Bau18] and others. 
A complete overview of nominal values for all dimensional parameters can be found 
in Table 6.2a, with exception of the rotational speed that is later varied systematically. 
It should of course be noted that, for instance, the nominal mass eccentricity is only 
relevant if an unbalance is considered at all and it can simply be ignored otherwise. 
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Table 6.1: Nominal parameter values used as reference for numerical analyses. 


(a) Dimensional formulation (available precision varies between parameters). 


Description Symbol Value Unit 
Bearing length L 38.10 mm 

3 Journal radius R 19.05 mm 

S Radial clearance C 50 pm 

Q Number of pads Np 3 = 
Gap orientation Xi We = 

~ Ambient pressure Poo 1013.25 hPa 

= Ambient temperature Doo 293.15 K 

™ Heat transfer coefficient Bo 100 W/(m?-k) 
Total number of bumps NpNg 24 — 

® Bump half-length bg 178 mm 

B Bump height hg 0.51 mm 

E Bump mass mg O1 g 

oi Bump stiffness kg 25 N/pm 
Friction coefficient uB 0.25 — 
Rotor mass MR 1 kg 
Mass distribution ER O1 — 

5 Mass eccentricity eR 7.5 pm 

2 Rotor shaft stiffness kr 10 N/pm 
Rotor damping coefficient dr 03 N-ms/pm 
Gravitational acceleration g 9.81 m/s? 


(b) Nondimensional formulation (rounded to three significant digits at most). 


Description Symbol Value Unit 
E Geometry number Ko Yao — 
Ö Number of pads Np 3 — 
Gap orientation x T2 = 
~ Ambient pressure (nondim.) Ps | 
= Ambient temperature (nondim.) To de = 
™ Thermal management number Qo 0.578 — 
Total number of bumps NpNg 24 — 
® Bump half-length (nondim.) Bg 35.6 — 
B Bump height (nondim.) Hpg 10.2 — 
E Bump mass number Mg 6.82 x 103 — 
N Bump frequency number Og 158 — 
Friction coefficient HB 0.25 — 
Rotor mass number MR 68.2 — 
Mass distribution ER 01 = 
5 Mass eccentricity (nondim.) ER 0.15 — 
2 Rotor frequency number OR 0.316 — 
Rotor damping number Dr 1.50 x 1072 — 
Gravitational loading number G 1.96 x 10 $ — 
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Before the dimensional parameter values from Table 6.2a can actually be inserted into 
Eqs. (5.4), (5-9), (5-10), (5.13), (5.15), and (5.18) to perform the nondimensionalization, 
an evaluation of Eq. (3.7) must yield the dependent parameter value u° ~ 11.6 Pa - s. 
Together with all the naturally dimensionless parameters, this finally results in the 
nominal set of nondimensional parameter values that Table 6.2b lists in its entirety. 
Reversely, the rotational speed being only expressed in the nondimensional form Ao, 
it reveals purposeful to explicitly calculate! the characteristic time scale tọ ~ 0.1 ms 
from Eq. (5.8), thus allowing to perform the back-transformation ng = Ao/(2ñto). 
As a rule of thumb for the interpretation of all following nondimensional results, 
this implies that a bearing number of Ag = 0.3 roughly corresponds to 30000 rpm. 


6.2 Application of Dynamical Systems Theory 


In the following introduction to mathematical fundamentals of dynamical systems, 
it is intended to always keep a close reference to the actual GFB-rotor model and to 
immediately apply the theory to problems with mostly nominal parameter values. 
Even though many effects could also be shown with simpler modeling approaches, 
all of the results presented hereafter are obtained assuming non-ideal gas behavior 
under non-isothermal conditions and taking into account foil compliance and friction. 


6.2.1 Classification of Solution Types 


For any possible starting point Sg € Q in the phase space of the dynamical system, 
which is specified as initial condition S(t = 1) = So at an arbitrary initial time Tọ, 
the solution T ++ S(t) of Eq. (5.27) directly describes the arising transient process. 
Nevertheless, the primary interest in theoretical and practical rotor dynamics often 


consists only in the long-term system behavior that can be observed for T — +00. 
In this context, the relevant analyses are usually concerned with limit sets toward 
which dynamical systems may evolve as time goes to positive or negative infinity. 
Topologically speaking, such limit sets Q+% C Q are subsets of the phase space 


and can be identified as points, curves, manifolds, or non-trivial fractal structures. 
More specifically, a limit set is called attractor (24 if it is reached asymptotically 
for T — + after starting from a surrounding region known as basin of attraction, 
while it is called repeller Q_. if it has that property for T + —oo (see, e.g., [MV11]). 


1 The found numerical value for the characteristic time scale remains valid throughout this entire chapter 
since p° and p° are basically constants for a given fluid and neither R nor C are to be varied anywhere. 
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2. 
(a) Fixed-point solution. (b) Periodic solution. (c) Quasi-periodic solution. 


Figure 6.1: Plots of transient rotor disk trajectories, which are obtained by numerical time integration and 
give a first impression of solution types toward which the dynamical system may evolve in the long term. 


In a corresponding numerical experiment, the first 200 ms of transient rotor drops are 
simulated by numerically integrating Eq. (5.27) over the time interval tT € [0,2000] 
using a fourth-order Rosenbrock method for systems of stiff ODEs (see, e.g., [HW96]). 
To create a simple starting point So, the fluid is initialized with ambient properties, 
while the foil structure and the centrally positioned rotor are undeformed and at rest. 
For the results, the dimension of the phase space being far too large to be visualized, 
the plots in Fig. 6.1 show projections? of the disk displacements XR(T) and Yr(T), 
which may evolve in completely different ways under slightly varied circumstances. 


Beginning with a rotational speed of Ap = 0.2 and setting the unbalance to er = 0 
as shown in Fig. 6.1a, the dynamical system and therefore the plotted displacements 
exhibit a transient evolution (dashed in green) that tends asymptotically toward a 
constant final state (green dot) and thereby reaches an attractive fixed-point solution. 
In Fig. 6.1b, when increasing the speed to Ag = 0.3 while still ignoring the unbalance, 
the transient motion (dashed in orange) now guides the dynamical system toward an 
attractive periodic solution or limit cycle (in orange), ending in self-excited vibrations. 
Finally, the addition of a mass eccentricity er = 0.15 with unaltered speed in Fig. 6.1c 
results in a quasi-periodic motion pattern (in blue), which is due to the occurrence of 
two incommensurable fundamental frequencies from self-excitation and unbalance. 
In fact, the presence of an unbalance does not necessarily induce quasi-periodicity, 
but can also lead to simply-periodic, multiply-periodic, or chaotic behavior [Boy11]. 
Since the present thesis is mainly concerned with fixed-point and periodic solutions, 
the reader is referred to the thesis of Becker [Bec20] on the topic of quasi-periodicity. 


? Alternative projections of other phase variables including densities, temperatures, displacements, and 
velocities could of course be plotted in direct analogy and would obviously lead to the same conclusions. 
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6.2.2 Stability and Bifurcation Analysis 


By pursuing the approach of brute-force time integration as in the previous section, 
it would be possible to get a more or less random overview of the solution behavior 
by just varying the rotational speed in small steps from one simulation to another. 
Regrettably, such a direct strategy turns out to be not only computationally inefficient, 
but principally unsuitable for several reasons when dealing with dynamical systems. 
First of all, time simulations can only give evidence for attracting limit sets (24, 
even though this shortcoming may be cured by additionally inverting the time axis T. 
Most importantly, however, it is commonly observed in nonlinear dynamical systems 
that there are two or more coexisting limit sets, each of which can then only be found 
by time integration when starting by chance within the associated basin of attraction. 
In this perspective, more systematic investigations may be conducted by making use 
of numerical continuation methods, which are essentially based on the observation 
that the continuous variation of a parameter constructs a family of related solutions. 
The following gives a brief explanation of the underlying theory and methodology 
based on some preliminary numerical results and without claiming to be exhaustive. 


For the initialization of a specific solution branch, a small initial bearing number ART 
is chosen such that there is a fixed-point solution SF! of Eq. (5.27) to fulfill (with i = 1) 


R{S", Af } = 0. (6.1) 


This is achieved either by trying if time integration runs into a non-evolving state or, 
more targeted, using an under-relaxed3 Newton-Raphson method (see, e.g., [FPS20]). 
From a practical point of view, finding a very first system state that solves Eq. (6.1) 
constitutes the most tricky part of each analysis since the Newton-Raphson method 
requires a good initial guess, which is more difficult the smaller the bearing number. 
Once a solution is known together with the corresponding bearing number, however, 
a predictor-corrector method can deliver a sequence of further fixed-point solutions 


(er) (Ge Ge ome (62) 
IG 


the existence of which is assured by the implicit function theorem (see, e.g., [Kuzo4]). 
For additional hints on the underlying procedure of pseudo-arc-length continuation, 
the interested reader is also referred to the documentation of AUTO 2000 [Doe+02]. 


Especially for practical applications, not only the existence and identification but also 
the stability of solutions is an important concern, which decides on whether some 


3 In order to improve the chance of convergence while not unnecessarily slowing down the calculation, 
the initial relaxation factor of 0.1 is increased successively during the final iterations until 1.0 is reached. 
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predicted behavior can actually be observed and is not of purely academic interest. 
In all of the following discussions, stability is understood in the sense of Lyapunov 
and therefore refers to perturbations in the initial conditions of a known solution. 
Denoting a perturbed solution by S,(T), this leads to the mathematical definition? 
that the solution of interest S(T) is called Lyapunov-stable if (see, e.g., [Hag78]) 


Ve>0, 36()>0: ||S.(t = T) — S(t = )]| < 6(e) 
= |S.(t)-S(t)||<e Vr>w. (63) 


For any arbitrarily small neighborhood> e > 0 around the solution S(T), this requires 
that all trajectories of perturbed solutions S,(T) that start at T = 1 from a sufficiently 
small neighborhood 6(¢) > 0 remain within the e-neighborhood for all times T > Tọ. 
For practical purposes, this definition of stability needs to be accompanied by some 
concept of attraction, which is introduced mathematically by (see, e.g., [Hag78]) 


3>0: [S(t =) — S(t =H) || <8 = lim ||S.() — $(x)|] =0. (64) 


Any solution that satisfies both Eq. (6.3) and Eq. (6.4) is called asymptotically stable 
and thereby represents an attractor (24% of the dynamical system (see, e.g., [NBg5]). 
Completing these definitions, a solution is called unstable if it is not Lyapunov-stable. 


In order to evaluate Eqs. (6.3) and (6.4) for the fixed-point solutions SF from Eq. (6.2), 
it is permissible by virtue of the Hartman-Grobman theorem to linearize Eq. (5.27) 
by simply calculating the Jacobians JẸ for all hyperbolic equilibria (see, e.g., [MV11]). 
As an efficient criterion, denoting by Af! the N complex eigenvalues of JR such that 


det (Jk = ART) Se sa Ne (6.5) 
it comes that a fixed-point solution SF is asymptotically stable if (see, e.g., [MV11]) 


RAR) < 0. 6.6 
Ps AU (6.6) 
In the contrary case, the solution of interest is principally assumed to be unstable 
even though there are undecidable cases in which the largest real part is just zero. 


A convenient way to make the behavior of a dynamical system visually accessible 
consists in the establishment of bifurcation diagrams, in which the progressions of 
selected solution-dependent quantities are plotted against a bifurcation parameter. 
With regard to rotating machinery, the displacements of the rotor are certainly among 


4 In definitions, some quantities such as £ and 6 are introduced independently of the actual nomenclature. 


5 The notion of neighborhood is defined in terms of a suitable norm, which may be the Euclidean norm. 
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the most natural variables to be observed when assessing the overall performance. 
To this effect, the bifurcation diagram in the upper part of Fig. 6.2a traces the extrema 
min; Yr(T) and max; YR(T), which gives the envelope of vertical rotor disk motions, 
as curves that are locally parameterized by the nondimensional rotational speed Ao. 
For any of the fixed-point solutions in Eq. (6.2), the disk remains at a specific place 
for all times and the represented envelope therefore degenerates into a single curve. 
In this particular plot, the nominal configuration according to Table 6.1 (in orange) 
is opposed to cases without (in green) and with increased rotor damping (in blue), 
which leaves the fixed-point solution branch in the middle completely unaffected®. 
As the bearing number grows, the rotor approaches asymptotically some position 
near the bearing center and thereby defines characteristic locus curves in Fig. 6.2b 
that are highlighted (in yellow) for the disk (top row) and the journals (bottom row). 
Moreover, all the plots in Fig. 6.2 contain information about the stability of solutions, 
which means that asymptotically stable fixed-point solutions according to Eq. (6.6) 
correspond to solid lines, whereas unstable solutions are shown as dashed lines. 
Obviously, there seems to be a critical bearing number AË for fixed-point solutions 
to become unstable, which is shifted to higher speeds the more the rotor is damped. 


When following the complex eigenvalues of the Jacobian in a closer vicinity of AHB, 


it is found that a pair of complex conjugate eigenvalues crosses the imaginary axis, 
which is a characteristic finding during the occurrence of a Hopf bifurcation (HB). 
Computationally, such HB points can efficiently be detected by introducing specific 
test functions to be evaluated as part of numerical continuation (see, e.g., [Doe+o2]). 
In addition to the loss of stability for fixed-point solutions, a HB also marks the birth 
of a branch of periodic solutions SP! (r), all of which satisfy Eq. (5.27) according to 


Egri (T) = RS" (r), AB}, ln) = S(r + Pi). (6.7) 
T 

For periodic solutions of autonomous dynamical systems, the period t” in Eq. (6.7) 
is a priori unknown and must thus be found by an appropriate solution strategy. 
Detailed descriptions of the orthogonal collocation method used here can be found 
in the book by Kuznetsov [Kuzo4], which also presents alternative shooting methods. 
Having calculated for Af! ~ A$P the first periodic solution SP! (t) with period t”!, 
a predictor—corrector method similar to that from Eq. (6.2) can provide a sequence 


(CU 7) ) _ (Geran), (S?2(x), AR, ms) (6.8) 


which is constructed by numerical continuation until a defined end point is reached 
or until convergence of the orthogonal collocation method fails (see, e.g., [Doe+o2]). 


iEN 


6 For better distinction, the three fixed-point solution branches are slightly shifted relative to each other. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.2: Results of numerical continuation and bifurcation analysis for different damping coefficients 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations/bifurcating solutions). 
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The assessment of stability according to Eqs. (6.3) and (6.4) is slightly more difficult 
for the periodic solutions SP! (t) from Eq. (6.8) than it is for fixed-point solutions. 
Using again the idea of linearization to judge the evolution of an initial perturbation, 
the Jacobian J}! (T) comes out as a time-dependent matrix with periodic coefficients. 
For dealing with the systems of linear ODEs with periodic coefficients that arise, 
the mathematical foundation is given by classical Floquet theory (see, e.g., [Hag78]). 
When calculating the fundamental matrix solution Mf’ (T) that is related to JẸ} (T) by 


MED = ROMEO), MRC = 0) = 1 (6.9) 
an evaluation after just one period yields the monodromy matrix Młi(t = T™'). 
The N complex eigenvalues AP of the monodromy matrix, which are also referred to 


as the Floquet multipliers of the periodic solution, can ultimately be determined by 
det [MR (T = T’) — ap] =0, n=1,...,N. (6.10) 


As a characteristic of autonomous dynamical systems, there is always one special 

Floquet multiplier Ae = 1 that represents tangential perturbations along an orbit, 

which can be visualized graphically by means of a corresponding Poincaré map. 

Based on the remaining (N — 1) Floquet multipliers, the Andronov—Witt theorem 
states that a periodic solution SP(T) is asymptotically stable if? (see, e.g., [MV11]) 

Pi 

en 


<1. (6.11) 


Similar to Eq. (6.6), the periodic solution is considered as being unstable otherwise, 
knowing that there may be again cases that cannot be decided by a linearized theory. 


When finally coming back to the bifurcation diagram in the upper part of Fig. 6.2a, 
upper and lower envelopes of vertical disk motions are obviously no longer identical 
for the periodic solutions that branch from the respective HB points (circle symbols). 
Also, the exemplary trajectories of the disk (top row) and the journals (bottom row) 
plotted in Fig. 6.2b now correspond to closed curves around each of the HB points, 
which illustrate the orbit® shapes in relation to the nominal clearance (dotted line). 
Due to the compliance of the foil structure, it is clear that the journals can travel to 
eccentricities e(t) > 1, which are even amplified for disk motions by the elastic shaft. 
By convention, asymptotic stability according to Eq. (6.11) is indicated by different 
line types for stable (solid lines) and for unstable (dashed lines) periodic solutions. 


7 Here, it is pointed out that |... | denotes the absolute value Y(R?(...) +3?(...)) of a complex number. 
8 Strictly speaking, the two-dimensional closed curves are not the actual orbits but only their projections. 
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As an additional information, the bifurcation diagram in the lower part of Fig. 6.2a 
visualizes all corresponding periods as nondimensional subsynchronous frequencies 


Ai = 27 
sub — Tri 


(6.12) 
in relation to Ag. This implies An /Ao & 0.5 for typical half-frequency whirling, 
whereas the frequency ratio would become AL / Ag = 1 for synchronous vibrations. 


Most remarkably, the variation of rotor damping has not only a quantitative effect 
on the amplitudes, but it also shows a qualitative influence on the system behavior. 
Beginning with the largest considered damping (in blue), a branch of stable periodic 
solutions emerges directly from the HB, which is therefore classified as supercritical. 
Practically speaking, the system may reach either a unique fixed-point solution or a 
unique periodic solution, depending on whether Ap is smaller or larger than ANB. 
For the nominal damping (in orange), however, the HB is found to be subcritical and 
thus yields a branch of unstable periodic solutions, along which Ao decreases again. 
Ultimately, it is only upon arrival at another critical bearing number ALP“ < AB 
that the occurrence of a limit point of cycles (LPC) results in a change of stability, 
which is due to a non-trivial Floquet multiplier crossing the complex unit circle at +1. 
At such LPC points (diamond symbols), the direction of continuation is reversed, 
which in the present case makes Ag grow until convergence fails (star symbols). 
Transferring this finding to applications, it is concluded that there are actually ranges 
of rotational speeds ALP“ < Ag < AH® in which stable fixed-point solutions and 
stable periodic solutions coexist. In such cases, the history of the system decides 
on the solution to be reached and large perturbations into the basin of attraction of 
another solution can cause sudden and radical changes during steady-state operation. 
A conservative design of rotating machinery thus requires a limitation of permissible 
rotational speeds to ALP“ instead of A}, which represents an important information 
that could impossibly be obtained by a purely linear GFB-rotor system analysis. 
Finally, things are even more complicated in the absence of damping (green lines), 
for which the HB is supercritical but this time followed by two subsequent LPCs. 
Looking closely at the bifurcation diagrams, there is not only a small region in which 
a fixed-point solution and a periodic solution coexist, but also a larger region with 
two coexisting periodic solutions of distinct vibration amplitudes and frequencies. 


With regard to the observed subsynchronous frequencies, the overall tendency shows 
that the fastest oscillations always occur near the place of birth of a stable limit cycle, 
whereas these vibrations slow down successively with increasing rotational speed. 
Nevertheless, the absolute level of frequency ratios is actually determined by the 
amount of damping and reaches from almost half-frequency whirling if undamped 
to approximately Vs-frequency whirling in the considered strongly damped scenario. 
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Transferring these findings to vapor-compression refrigeration machinery on GFBs, 
it should be remembered that the rotor damping number DR is meant to represent 
viscous forces that the surrounding working fluid exerts on a laterally moving rotor. 
Since vapor-liquid phase transitions are usually accompanied by an important rise 
in viscosity as already shown in Fig. 3.2e, the crucial influence of this parameter 
reveals in fact how sensitive such systems are to varying thermodynamic conditions. 
In other words, the overall machine seems much more prone to self-excited vibrations 
when the rotor is immersed in pure vapor rather than in a vapor-liquid mixture, 
but such a conclusion cannot be made before studying consequences for the GFBs. 


6.3 Particular Features of Refrigerant Lubrication 


It is not surprising that the essential difference between refrigerant-lubricated GFBs 
and air-lubricated GFBs consists in the possible occurrence of phase transitions. 
However, even though this effect is included in the comprehensive fluid modeling, 
none of the results presented so far allow for an actual discussion of consequences. 
As a general rule, condensation becomes an issue when the fluid is compressed to 
the saturated vapor pressure, which is reached the earlier the lower the temperature. 
Intuitively, it can be expected that this often leads to scenarios in which only the 
heating during operation prevents the vapor from condensing under heavy loading. 
In this context, the thermal management number ©) is a decisive design parameter of 
the GFB model because it describes the ability of transferring heat to the environment. 
Concerning the characterization of heavy loading, there should be a clear distinction 
between steady-state operation, in which the weight of the rotor is most important, 
and situations in which dynamic forces due to self-excited vibrations are dominant. 
In the following part of the analysis, fixed-point and periodic solutions are therefore 
investigated with a focus on thermodynamic properties that describe the refrigerant. 


6.3.1 Fluid Behavior and Phase Transitions 


Looking first at systems that operate in an equilibrium position, the expected effects 
can best be emphasized by considering a rather heavy rotor with MR = 341 (5kg). 
Leaving the other parameters from Table 6.1 unchanged, which implies @g = 0.578, 
this set-up is opposed to one with more effective cooling behavior such that @5 = 144. 
At the rotational speed Ag = 0.5, the dynamical system for the heavy rotor possesses 
one stable fixed-point solution for either of these thermal management numbers. 
Since the discretized densities and temperatures are part of any solution state vector, 
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from which other thermodynamic properties can be retrieved via equations of state, 
it is straightforward to plot the corresponding fields in Fig. 6.3 through interpolation. 


The temperature field for the nominal thermal management in Fig. 6.3e shows that 
the lubricant is subject to a noticeable heating from the ambient temperature of 20 °C 
to approximately 60°C along the pad that is closest to the displaced shaft journal. 
Also, the temperature decreases just very slightly toward the axial bearing edges, 
which means that the saturated vapor pressure is highly above the ambient pressure 
throughout practically the entire active zone? of the bearing. The discontinuities at 
all fixation gaps reflect the simplifying assumption of inflow and outflow boundaries 
without mixing, however, this does not much affect temperatures in the active zone". 
Considering the corresponding density field plotted in Fig. 6.3c, the most important 
compression occurs in the converging region where also the temperature rise begins, 
whereas a drop to subambient densities is found toward the trailing edge of that pad. 
As a matter of fact, the observed expansion partly explains why the hot refrigerant 
transported through the active zone suddenly stagnates at a certain temperature level. 
Finally, the calculated pressure field according to Fig. 6.3a clearly reveals that the 
peak pressure reaches not even half the locally prevailing saturated vapor pressure, 
meaning that no vapor-liquid phase transition is predicted for the considered case. 
Obviously, since the pressure satisfies Dirichlet boundary conditions at all pad edges, 
it may be represented as one continuous field that is not differentiable at the gaps. 


As can be seen from the temperature field for better thermal management in Fig. 6.3f, 
the lubricant in the active zone is heated not much above the ambient temperature. 
Most importantly, this implies that the saturated vapor pressure is not much higher 
than ambient pressure either, thus making the refrigerant prone to phase transitions. 
Indeed, the density field in Fig. 6.3d shows the fluid compressed to an extent that 
the corresponding pressure field in Fig. 6.3b might seem disproportionately small, 
however, the overall pressure magnitude is of course prescribed by the rotor weight. 
Rather, the presence of vapor—liquid mixture locally prevents any pressurization 
beyond the saturated vapor pressure, which effectively softens the lubricant film. 
Since the limiting saturated vapor pressure actually varies only with temperature, 
this straightly cuts off pressure peaks under perfectly isothermal conditions [LBS17c]. 
Finally, looking closely at both the pressure distribution and the density distribution, 
the two-phase flow collapses at a certain place more abruptly than when it is created. 
This effect weakens with increasing regularization of the latent heat of condensation 
as it is mostly due to discontinuities in heat capacity and compressibility coefficient. 


9 Active zone refers to that region of the lubrication gap in which most of the pressure build-up occurs. 
1° This anticipates what is found in Fig. 6.5e for single-pad GFBs without fixation gap in the active zone. 
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Figure 6.3: Plots of thermodynamic properties observed during steady-state operation of three-pad GFBs 
with mediocre thermal management (left column) or with more effective cooling behavior (right column). 
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Figure 6.4: Plots visualizing condensation during steady-state operation of three-pad or single-pad GFBs. 


When comparing the two pressure distributions in Figs. 6.3a and 6.3b to each other, 
they differ in shape but lead to identical bearing forces that support the rotor shaft. 
As the most important difference, phase transitions flatten the pressure field and 
thereby expand the active zone toward the axial bearing edges, which is accompanied 
by an increased journal eccentricity of € ~ 0.95 versus € ~ 0.85 without condensation. 
In one way or another, this observation demonstrates how phase changes may have 
a detrimental influence on the steady-state performance of GFB-rotor systems. 
Firstly, the finding of systematically increased eccentricities means that the minimum 
permissible film thickness condition according to Eq. (2.17) is violated more easily, 
which signalizes a potentially reduced load-carrying capacity of concerned bearings. 
Secondly, there is no guarantee that limited pressure magnitudes due to condensation 
can actually be compensated at all by pressure redistribution and journal relocation. 
With the vapor quality field in Fig. 6.4a, which belongs to the right column of Fig. 6.3, 
it becomes evident that the considered steady-state operating conditions represent 
a case with the compensation capabilities of the system being practically exceeded. 
More precisely, pure vapor is only present in regions that are largely insignificant for 
the performance and two-phase flow is predicted for almost the entire active zone. 


For the sake of completeness, a similar analysis is conducted for single-pad GFBs 
with the results in Fig. 6.5 plotted for the same two thermal management numbers 
and Fig. 6.4b showing the vapor quality field for the scenario with improved cooling. 
In order to not repeat all the observations that are comparable to three-pad GFBs, 
only the most striking differences and some additional details are addressed hereafter. 
Most importantly, given that there is only one fixation gap located at ọ = x; = 71, 
the heated lubricant is transported farther away from the active zone while basically 
maintaining constant temperatures near the bearing centerline as Fig. 6.5e suggests. 
It can be concluded that fluid exchange with the environment via the fixation gaps 
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Figure 6.5: Plots of thermodynamic properties observed during steady-state operation of single-pad GFBs 
with mediocre thermal management (left column) or with more effective cooling behavior (right column). 
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has a surprisingly little effect on the overall temperature level in the active zone. 
Anyway, the reduced number of fixation gaps renders the bearing more flexible 
in the sense that it can adapt better to phase transitions by redistributing pressure 
without being constrained by pressure losses at the circumferential pad boundaries. 
Consequently, single-pad GFBs have more reserves to accommodate two-phase flow 
and generally exhibit smaller amounts of condensed liquid under similar conditions, 
which becomes obvious when the two vapor quality fields in Fig. 6.4b are opposed. 


Until here, the analyses of three-pad and single-pad GFBs in steady-state operation 
are all performed for just one exemplary rotational speed of Ag = 0.5 (50000 rpm), 
so that it remains to study the influence of the bearing number on the fluid behavior. 
As a reference for thermodynamic properties along the centerline of three-pad GFBs, 
the plots in Figs. 6.6a and 6.6b originate from the respective fields in Figs. 6.3 and 6.4, 
which are completed by the associated progressions of viscosity and vapor quality. 
Similarly, the stable fixed-point solutions that are found for Ag = 1.0 (100000 rpm) 
lead to the representations in Figs. 6.6c and 6.6d, while those in Figs. 6.6e and 6.6f 
correspond to Ag = 1.5 (150 000 rpm) with different thermal management numbers. 
In all of these plots, the pressure distribution as predicted by the perfect gas model 
from Section 3.1.3 is also indicated to subsequently allow for direct comparisons. 


Looking at Figs. 6.6a, 6.6c, and 6.6e to discuss the nominal case where @ = 0.578, 
the variation of rotational speed most obviously affects temperatures (in magenta). 
Since the pad in the active zone features no more region of subambient pressures 
for higher rotational speeds, the heating of the fluid is no longer counteracted by 
expansion and temperatures grow further until the flow arrives at the fixation gap. 
For the distant pads, the influence of rotational speed becomes even more significant, 
but the overall performance is not much concerned by temperatures in those regions. 
Most importantly, no phase transitions are predicted for any of the considered Ag, 
resulting in virtually no difference between the pressure profiles that are obtained 
from the real gas model (solid in green) and the perfect gas model (dashed in green). 
Furthermore, the nondimensional viscosity (in yellow) coincides almost perfectly** 
with the nondimensional temperature for the trivial vapor quality W = 1 (in black). 


For the more effective cooling with Oo = 144 as shown in Figs. 6.6b, 6.6d, and 6.6f, 
there is an interesting relationship between rotational speed and two-phase flow. 
In fact, the amount of condensed liquid decreases with increasing bearing number, 
which cannot be explained by the observation of just slightly higher temperatures. 
Instead, it is the characteristic shape of the pressure distribution that exhibits a 
pronounced peak for slower speeds and becomes much wider for higher speeds 


™ For better distinction, temperature curves and viscosity curves are slightly shifted relative to each other. 
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Figure 6.6: Plots of thermodynamic properties observed during steady-state operation at different speeds 
with mediocre thermal management (left column) or with more effective cooling behavior (right column). 
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when the shaft journal approaches the bearing center as already seen in Section 6.2.2. 
Moreover, phase transitions can be accompanied by tremendous increases in viscosity, 
which is very sensitive to the vapor quality but partially compensates pressure losses. 
Finally, the pressure profiles in Figs. 6.6b and 6.6d clearly illustrate the limitations of 
the perfect gas model, which is not capable of reproducing phase transitions and 
may therefore drastically overestimate peak pressures and load-carrying capacities. 


In conclusion, vapor-liquid phase transitions of the lubricant become more critical 
the smaller the rotational speed in steady-state operation of refrigeration machinery. 
This seems to be an important complement to the original findings of Garcia [Gar12], 
who shows that the amount of condensed liquid increases with the rotational speed 
if a shaft journal is kept fixed at constant eccentricity under stationary conditions. 
Here, taking into account the interaction with a rotor, the trend is just the opposite 
since the reduced eccentricities at higher speeds do not cause sharp pressure peaks. 
Apart from this, the found results reveal a remarkable robustness and adaptability 
of refrigerant-lubricated GFB-rotor systems, making it rather difficult to identify 
scenarios in which condensation is not entirely prevented by increasing temperatures. 
Therefore, knowing that vapor-liquid phase transitions are usually not beneficial, 
overly effective cooling of refrigerant-lubricated GFBs may turn into a disadvantage. 
Also, care must be taken in the design of such machines that inflowing refrigerant 
from the bearings’ environment is not too cold or even partially condensed already. 


6.3.2 Significance to Dynamic Performance 


In most situations, the critical bearing number AHP for the birth of periodic solutions 
falls into the region of high rotational speeds in which only pure vapor is present. 
However, this says nothing about the further evolution of fluid behavior in particular 
along branches of periodic solutions that finally lead to large vibration amplitudes. 
Since all discretized densities and temperatures are part of the overall state vector 
just like the rotor displacements, it does not surprise that for any periodic solution 
the closed curve of the orbit can simply be projected on combinations of D; j and T; j. 
When interested in vapor-liquid phase transitions, this bears the practical advantage 
that the saturated vapor line is also defined in terms of density and temperature, 
which divides such subspaces of Q into a vapor (V) and a vapor-liquid (V-L) region. 
For the plots in Fig. 6.7, the fluid is observed at three locations on the active pad: 
p = ~"/12 (near leading edge), p = 7/6 (center), and ọ = 57/12 (near trailing edge). 
Again, the left column corresponds to the nominal thermal management number, 
whereas the right column is now obtained for ©ọ = 12.4, which is deemed sufficient 
to represent effective cooling for the lightweight nominal rotor with Mr = 68.2 (1kg). 
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Figure 6.7: Plots of densities and temperatures observed during self-excited vibrations at different speeds 
with mediocre thermal management (left column) or with more effective cooling behavior (right column). 
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When performing whirling motion, the shaft journals circulate counterclockwise 
within the bearings while typically coming rather close to the bounding foil structure. 
Fixing an observer to an arbitrary circumferential position in the lubrication gap, 
subsequent passages of the journal are separated by the period of a periodic solution. 
Here, the density—temperature cycles in Fig. 6.7 are not passed through uniformly 
because most of the time the fluid remains close to ambient conditions (black dot) 
and then experiences only one short moment per period of heavy dynamic loading. 
In this connection, it is observed that the approaching journal is always preceded by 
a fluid depression zone with subambient densities but superambient temperatures, 
which becomes obvious when following the projected orbits in clockwise direction. 


Beginning with nominal thermal management, the largest fluctuations in densities 
and temperatures are found in Figs. 6.7c and 6.7e, depending on the rotational speed. 
For the largest considered bearing number of Ag = 0.5 (50 000 rpm), this implies that 
the temperature oscillates between 20°C and approximately 140 °C in that region, 
which is mostly due to the transport of heated refrigerant along the lubrication gap. 
Similar to the situation in steady-state operation, this shows once more how any 
compression of the fluid is always accompanied by high temperatures so that the 
thermodynamic state nowhere ever crosses the saturated vapor line (dashed in gray). 
When lowering the level of reached temperatures through a more effective cooling, 
all of the significantly flattened cycles in Figs. 6.7d and 6.7f clearly include transitions 
into the two-phase region without a quantifiable dependency on the rotational speed. 
Altogether, this reveals a second regime in which condensation plays a crucial role, 
this time restricting the speeds for which self-excited vibrations can be sustained. 


To further clarify the relevance of this effect, numerical continuation is performed 
with different thermal management numbers, for which Fig. 6.8 shows the outcome. 
Principally, there seems to be no important discrepancy between the established 
bifurcation diagrams in Fig. 6.8a, in particular as long as only pure vapor is present. 
Along the fixed-point solution branches, the impact of varying cooling performance 
on rotor shaft displacements is so marginal that it remains virtually invisible here”. 
With regard to periodic solution branches, however, the effect of pressure limitations 
through condensation is manifested by a slight increase in vibration amplitudes. 
Moreover, believing numerical convergence failure to correlate with a collapsing film, 
there is a clear tendency toward solid-to-solid contact the more condensation occurs. 
Apart from this, however, there is not much of a difference between the resulting 
rotor trajectories that are represented in Fig. 6.8b for the sake of completeness only. 


12 The classification into clockwise and counterclockwise motion refers to the representations in Fig. 6.2b. 
3 This explains why the plots in Section 6.3.1 require a heavier rotor than the nominal one to be illustrative. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.8: Results of numerical continuation and bifurcation analysis for different thermal management 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations /bifurcating solutions). 
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6.4 Influence of Foil Structure Design Parameters 


Using the established model as a virtual prototype for general GFB-rotor systems, 
an overwhelming variety of numerical parameter or sensitivity studies is imaginable. 
For the present thesis, which focuses here on the qualitative rotordynamic behavior, 
the selected design parameters are those concerned with the pad configuration and 
with properties of the foil structure, while further studies are published elsewhere”. 


6.4.1 Number and Orientation of Pads 


Prior to dealing with foil properties in the narrower sense, it is worthwhile to study 
how the global pad configuration itself affects the dynamic behavior of the system. 


According to Table 6.1, the nominal three-pad GFBs are oriented by x1 = 7/2 (90°), 
which is exactly the angle that is also used in Fig. 2.1 for the exemplary sketches. 
In the bifurcation diagrams in Fig. 6.9a, this gives just the same curves (in orange) 
that serve as a reference in Fig. 6.2a and show a subcritical HB followed by an LPC. 
Rotating the pads by small steps of 30° to four distinguishable constellations in total, 
only x1 = 27/3 (120°) yields similar curves (in blue) to those of the nominal setting. 
Notably, the reorientation of fixation gaps lowers the critical bearing number APE 
by approximately 0.07 (7000 rpm), whereas the detected HB remains subcritical. 
Anyway, since the arising branch of unstable periodic solutions is now stabilized 
earlier by an LPC, the critical bearing number ADP Cis finally similar for 120° and 90°. 
Obviously, things are getting much more interesting in situations where one of the 
fixation gaps interferes with the active zone of the bearing for a given load direction. 
Since the escaping lubricant hinders local pressure build-up, steady-state operation 
is only possible when the attitude angle and the eccentricity adapt accordingly. 
For x1 = °”Ys (150°), no such position at all is found for extremely low speeds and 
the fixed-point solution branch (in magenta) begins only after the shaft journal spins 
fast enough to eventually climb over the gap. Also, the critical bearing number AFP 
is reduced further by another 0.04 (4000 rpm) to encounter a subcritical HB, however, 
there is not much overlap where stable fixed-point and periodic solutions coexist. 
Along the branch of stable periodic solutions, amplitudes initially remain small but 
grow abruptly at bearing numbers for which the nominal system just starts to vibrate. 
Coming finally to the seemingly most extreme gap orientation angle x; = 7/3 (60°), 
there are two branches (in green) of partially coexisting stable fixed-point solutions, 
interrupted by an unstable branch due to the double occurrence of a limit point (LP). 


™4 A comprehensive list of publications related to this thesis can be found directly after the bibliography. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.9: Results of numerical continuation and bifurcation analysis for different gap orientation angles 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations /bifurcating solutions). 
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Later, the HB that is reached at approximately half the critical bearing number AFP 
of the nominal case turns out to be supercritical and directly produces a branch of 
stable periodic solutions. Once more, the amplitudes of these self-excited vibrations 
remain rather small and may be tolerable until finally approaching the other curves. 
From the lower part of Fig. 6.9a, it becomes obvious that such early vibrations occur 
at higher frequency ratios (in green and magenta) as later subsynchronous vibrations. 
Lastly, some exemplary fixed-point and periodic solutions are visualized in Fig. 6.9b 
to illustrate how rotor locus curves (in yellow) are relocated or cut into two parts, 
whereas the orbit shapes exhibit a remarkable sensitivity to the pad orientation. 
Looking in particular at displacements of the disk (top row), it seems like the rotor 
bounces from one pad’s center to the next, resulting in almost triangular trajectories. 


Going one step further, the number of pads Np represents certainly one of the most 
fundamental parameters for describing full-arc, single-pad, and multi-pad bearings. 
With regard to this investigation, it should be noted that the number of bumps Ng 
is always adjusted such that NpNg = 24 for the sake of comparability of results. 
Nevertheless, it is deemed reasonable to select the most natural orientation angle x1 
for each of the respective designs as indicated by the plots in Fig. 6.10b (bottom row). 
Moreover, since the bump foil model from Chapter 4 is not applicable’ for Np = 0, 
this particular case must fall back on the description Q(¢,T) = 0 of rigid bearings. 
Considering the bifurcation diagrams in Fig. 6.10a, there is an obvious relationship 
between the number of pads and the observed behavior of the GFB-rotor system. 
As a general tendency, the critical bearing number Ah® clearly increases with Np, 
however, the found HB changes from supercritical for full-arc or single-pad bearings 
to subcritical for multi-pad bearings. Consequently, the gain in stability may be 
reduced when aiming at a conservative machine design with ALC the speed limit. 
Comparing periodic solutions for self-excited vibrations, there is not much difference 
in amplitude but the frequency ratio is visibly lower the larger the number of pads. 
On the other hand, shaft journal eccentricities are much higher for multi-pad GFBs, 
which ultimately implies a reduced load-carrying capacity in steady-state operation. 
Coming to the representative trajectories of rotor disk and shaft journals in Fig. 6.10b, 
it is interesting to see how the foil fixation gaps are reflected in the orbit shapes and 
how the absence of such gaps prevents the journal from coming close to the top foil. 


Altogether, it can be concluded that the configuration of pads gives the developer 
decisive influence on the dynamic performance of refrigeration machinery on GFBs. 
Depending on the requirements, it might either be beneficial to expand the range for 
steady-state operation by using multi-pad GFBs or it might be preferable to ensure 


"5 In fact, this is due to the actual design of bump-type foils and not a fault of the modeling approach. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.10: Results of numerical continuation and bifurcation analysis for different numbers of pads 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations /bifurcating solutions). 
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Figure 6.11: Plots visualizing top foil contours during steady-state operation of three-pad GFBs. 


higher load-carrying capacities at the cost of less stability with single-pad GFBs'f. 
In either way, it is important to pay attention to the orientation of the foil structure 
relative to the predominant load direction, which quickly becomes a delicate task 
when facing multi-pad GFBs with foil fixation gaps around the entire circumference. 


6.4.2 Foil Properties and Dry Friction 


Among the few parameters to which the considered foil structure model is reduced, 
it is mainly the bump stiffness and the friction coefficient that affect the entire system. 
To get a principal idea of foil structure deformations during steady-state operation, 
which can be reconstructed from bump displacements in the solution state vector, 
the plots in Fig. 6.11 show stationary top foil contours’” for Ag = 0.1 (10 000 rpm). 
Comparing the results in Fig. 6.11a that are obtained for a rotor on rather soft GFBs 
to those for the nominal configuration in Fig. 6.11b, the plane deformation fields 
clearly differ not only in magnitude but also qualitatively along the circumference. 
Principally, this is due to the rotor adapting to increased foil structure compliance 
with higher shaft journal eccentricities accompanied by smaller attitude angles, 


which leads to larger regions of subambient pressures and positive deformations"®. 


For systematic numerical continuation, the frequency number Qg = 15.8 (2.5 N/um) 
of the nominal case is lowered to Qg = 10.0 (1.0 N/pm) and Qg = 7.1 (0.5N/pm), 
which leads to the bifurcation diagrams that Fig. 6.12a shows contrasted to Qg — oo. 


16 With the same benefits, there are also leaf-type GFBs [Agro7] that realize compliant full-arc bearings. 
17 The top foil fixations are only indicated for the purpose of illustration without being actual deformations. 
18 Even though not explicitly considered in the modeling, this can be interpreted as foil structure lift-off. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.12: Results of numerical continuation and bifurcation analysis for different bump frequencies 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations/bifurcating solutions). 
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While the critical bearing number A$P increases just slightly with decreasing Qp, 
other tendencies are quite obvious and reveal that the region of coexisting solutions 
due to a subcritical HB reacts very sensitively to the stiffness of the foil structure. 
From this observation, it can be concluded that structural compliance is generally 
beneficial for avoiding the sudden onset of self-excited vibrations during operation, 
but always accompanied by larger eccentricities and reduced load-carrying capacities. 
At the same time, the frequency ratio clearly decreases together with Qg and finally 
results in almost Y/s-frequency whirling for soft structures at high rotational speeds. 
Surprisingly, the rotor disk orbits (top row) in Fig. 6.12b are not that much affected 
even though the effect on shaft journal trajectories (bottom row) appears significant, 
so less deformation of the foil structure means more deformation of the rotor shaft. 


As indicated among the nominal parameter values, all previously presented results 
assume a friction coefficient of ug = 0.25 for the regularized Coulomb friction model. 
Since friction is generally regarded as being essential for the functioning of GFBs, 
the influence of this parameter on the system behavior merits a closer investigation. 
To this purpose, the bifurcation diagrams in Fig. 6.13a oppose friction coefficients 
ranging between ug = 0.15 and ug = 0.75 to the frictionless limiting case ug = 0.00. 
Obviously, fixed-point solutions remain completely unaffected as long as the friction 
coefficient is incorporated into the dynamical system via such a regularized model. 
On the other hand, considering a bristle model that truly accounts for static friction, 
the found ambiguity between multiple fixed-point solutions reveals largely irrelevant 
according to numerous time simulations’? performed with different initial conditions. 
Notwithstanding this finding, the question of how sticking contacts might affect ARE 
can hardly be answered with regularized friction forces that vanish at zero velocity. 
To roughly assess the matter, there is evidence that bumps that stick at both sides 
appear significantly stiffer toward the fluid than bumps that are free to slide [Alv18]. 
Relating therefore the influence of static friction to an increase of bump stiffness, 
it is already known that the stability of fixed-point solutions is insensitive to Qg. 
Consequently, the actual benefits of friction become only evident after vibrations 
emerge with amplitudes that are generally larger the smaller the friction coefficient, 
which can easily be seen in Fig. 6.13b from the respective orbits of periodic solutions. 
However, this reduction of vibration amplitudes remains somewhere limited and 
therefore cannot be amplified infinitely by further increasing the friction coefficient, 
so that ug = 1.00 (not depicted) would not differ much from ug = 0.75 (in magenta). 
On the contrary, the beneficial effect can even reverse if the coefficient is so high that 
it impedes relative motion in the contacts to an extent that less energy is dissipated, 
which makes amplitudes grow again as can be seen in Fig. 6.13a for Ag = 0.45. 


*9 It should be recalled from Section 4.4 that bristle models cause some difficulties for bifurcation analyses. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.13: Results of numerical continuation and bifurcation analysis for different friction coefficients 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations/bifurcating solutions). 
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Finally, it is noted that the subsynchronous frequencies of self-excited vibrations are 
basically independent of the friction coefficient, which corresponds pretty much to 
the expectations when adding Coulomb damping (not viscous damping) to a system. 


Now that the influence of frictional dissipation on the overall bearing-rotor system 
is known, a closer look at the friction forces themselves can give additional insights. 
To this purpose, the representative periodic solution found at Ag = 0.3 (30000 rpm) 
for the nominal system is calculated while tracing the progressions of friction forces, 
which is done during an arbitrary two-period?° time interval such that t € [0,2t°"]. 
Since there is a large number of NPNg = 24 contacts with as many friction forces n,m, 
the plots in Fig. 6.14 just consider the contacts?" m = 1,5,8 of the active pad n = 3. 
Performing the very same procedure with the elasto-plastic bristle friction model 
according to Dupont et al. [Dup+o2] that is discussed and described in Section 4.4, 
this reference model (in green) can be contrasted to the regularization (in magenta). 
Depending on the nature of friction modeling, predicted transitions (dashed in gray) 
between forward/backward sliding and sticking must be interpreted differently. 
While the smoothed sign function of the regularized friction model only indicates 
the direction of motion by +1 and tends to zero where its validity range is exhausted, 
a bristle switch function remains close to zero before and during contact breakaway. 
Concerning the bristle switch function, it may surprise that especially in Fig. 6.14b 
there seems to be no reverse motion, which is contradictory for periodic solutions. 
Hence, it must be clarified that backward sliding does occur when friction models 
are disregarded and friction forces are set to zero because of subambient pressures. 


In general, both models are in good agreement during phases of forward sliding 
although the regularization slightly overestimates the temporal extension of peaks. 
During backward sliding, the regularized model reveals less accurate, which can be 
explained by the fact that such motions mostly result from bump-bump interaction. 
In fact, when bump displacements propagate through a foil strip, sticking contacts 
have an isolating effect whereas non-sticking models overestimate the foil mobility, 
which becomes even worse in Figs. 6.14e and 6.14f where a contact is entirely stuck. 
Nevertheless, this is unlikely to alter the order of magnitude of dissipated energy, 
which after all is mainly determined by the level of pressure-induced normal forces. 
Altogether, the strong fluctuation of normal forces and therefore of friction plays an 
important role that is mostly ignored by studies using equivalent viscous damping. 


20 Here, it should be remembered that T”! from Eq. (6.7) denotes the period of a periodic solution SP” (T). 


21 Even though clear tendencies are generally observed from the center of a pad to the respective edges, 
there can also be particular cases with surprising discrepancies between neighboring contacts [LBS17a]. 


94 


6.4 Influence of Foil Structure Design Parameters 


Friction Force 773 1 


(a) Active pad’s leading edge (regularized friction). 


(e) Active pad’s trailing edge (regularized friction). 


Friction Force 713 5 


Friction Force 7738 


0.02 — mm u 
1 ft I 
Ioi It 
1 | bo 
GOT ries eve Lead tiie 
iN IN 
1] \ | | 
0.00 0 
i | ro 
1 | 1 | 
0.01 Hohe. EOE 
I I 
‘ae N! 1 
0.02 LU 5 À M I 4; 
0 50 100 150 


Nondimensional Time T 


0.02 +1 
0.01 
0.00 0 
ry 
—0.01 I i 
oie \ı 
11 | | 
0.02 og 
0 50 100 150 


Nondimensional Time T 


(FEx)uss pozuemnsoy 


(££y)u8s paziepnSay 


(c) Active pad’s center (regularized friction). 


+1 


0.02 me | 
I 1 \ 
1 | IA 
0.01 + Le | 4 
El IA 
i i 1 1 
0.0 u U | 0 
REP TY 
001 Len 
1 1 
LE U il 4 
0.02 ii L Ka i 
0 50 100 150 


Nondimensional Time T 


(S€y)u8s paziepnSay 


Friction Force 773,1 


rm 


uorpung MS dIs-PuS 


rm 


0.02 77 — + 
ji 11 
| | | | 
0.01 L . li i 
lal: B A 
0.00 = ER vn 0 
—0.01 + 
0.02 j i i 
0 50 100 150 


Nondimensional Time T 


(b) Active pad’s leading edge (bristle friction). 
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(d) Active pad’s center (bristle friction). 
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Figure 6.14: Plots of friction forces and stick-slip switch functions observed during self-excited vibrations 
assuming either regularized Coulomb friction (left column) or elasto-plastic bristle friction (right column). 
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6.5 Remarks on the Role of Rotor Unbalance 


Looking from the point of view of rotor dynamics at results presented and discussed 
throughout this chapter, probably the most questionable idealization results from the 
assumption of perfectly balanced rotors. In real-world rotating machinery, however, 
there is of course always some residual unbalance despite all technical progress, 
which makes it necessary to see how the previously drawn conclusions are affected. 
First of all, the presence of an unbalance implies that the dynamical system possesses 
no fixed-point solutions since there are no actual resting positions of fluid, structure, 
and rotor. Instead, it can be supposed that the system state fluctuates periodically 
around the former equilibrium point, depending on the amount of mass eccentricity. 
In numerical continuation, a common way to find such periodic solution branches are 
so-called homotopy methods that temporarily consider er the bifurcation parameter. 
Starting from a fixed-point solution of the corresponding balanced system (er = 0), 
this creates a continuous transformation to the respective er > 0 (see, e.g., [Doe+o2]) 
and thereby initializes a periodic solution branch to be continued by variation of Ag. 


Referring to the perfectly balanced nominal system (in orange), the plots in Fig. 6.15 
show the effect of different mass eccentricities by tracing the arising branches of 
periodic solutions (in green, blue, and magenta) around former fixed-point solutions. 
In analogy to scenarios without unbalance, these periodic solutions are only stable 
for small Ag and become unstable beyond AN®® according to the Floquet theory. 
Interestingly, there is a clear tendency of critical bearing numbers increasing with 
mass eccentricity, which in relation to A$P shows a stabilizing effect of unbalance 
at the cost of larger vibration amplitudes that finally cause failure (in magenta). 
At APE, two complex conjugate Floquet multipliers that cross the complex unit circle 
indicate the occurrence of a Neimark-Sacker bifurcation (NSB) or torus bifurcation, 
which is supposed to result in the birth of quasi-periodic solutions (see, e.g., [MV11]). 
Since numerical continuation along such branches is the subject of ongoing research 
and as it reveals delicate to apply according methods to large dynamical systems, 
the bifurcation diagrams in Fig. 6.15a and orbits in Fig. 6.15b ignore these solutions. 
Instead, time simulations”? are carried out for discrete bearing numbers Ag > Aa 
with a moderate mass eccentricity er = 0.15 (7.5 um) to then illustrate in Fig. 6.16 
the almost bizarre variety of attractors Q+% in this regime of the dynamical system. 


>? For the sake of better visibility, the displayed trajectories correspond to relatively short periods of time 
as otherwise densely filled areas without distinguishable lines would result for quasi-periodic solutions. 
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(a) Bifurcation diagrams for vertical rotor disk displacements (top) and frequency ratios (bottom). 


(b) Representative trajectories of the rotor disk (top row) and of the shaft journals (bottom row). 


Figure 6.15: Results of numerical continuation and bifurcation analysis for different mass eccentricities 
(solid lines: stable solutions, dashed lines: unstable solutions, markers: bifurcations/bifurcating solutions). 
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6 Nonlinear Effects in GFB—Rotor Systems 


(b) Rotational speed Ag = 0.325. (c) Rotational speed Ap = 0.350. 


(d) Rotational speed Ag = 0.375. (e) Rotational speed A5 = 0.400. (f) Rotational speed Ag = 0.425. 


(g) Rotational speed Ap = 0.450. (h) Rotational speed Ap = 0.475. (i) Rotational speed Ag = 0.500. 


Figure 6.16: Plots of transient rotor disk trajectories for some Ag > ANS with mass eccentricity er = 0.15. 
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7 Conclusion 


7.1 Summary 


The present thesis addresses the dynamics of rotors on refrigerant-lubricated GFBs 
and establishes a new modeling approach to enable stability and bifurcation analyses. 
Owing to the precisely identifiable fluid-structure-rotor interaction mechanisms, 
three submodels that are of reasonable complexity but nevertheless take into account 
all relevant nonlinearities can finally be transformed into a single dynamical system. 
As an essential model feature, the non-ideal characteristics of a typical refrigerant 
that may undergo vapor-liquid phase transitions are described by thermodynamic 
equations of state to be included in modified Reynolds and temperature equations. 
Also, it becomes feasible with the proposed lumped-element foil structure model 
to include dry friction in various ways reaching from highly efficient regularized 
Coulomb friction models to bristle models that can reproduce stick-slip transitions. 
Altogether, this makes the entire problem accessible to rigorous mathematical theory 
and allows for developing a monolithic research code with interchangeable modules. 


When studying the long-term behavior of refrigerant-lubricated GFB-rotor systems, 
transient states evolve to fixed-point, periodic, quasi-periodic, or chaotic solutions 
depending on design parameters and on the systematically varied rotational speed. 
For perfectly balanced rotors, it is found that equilibrium positions lose their stability 
via subcritical or supercritical Hopf bifurcations associated with some critical speed, 
which typically results in an immediate onset of undesirable self-excited vibrations. 
Under the arising dynamic loading, refrigerant-lubricated GFBs become particularly 
susceptible to phase transitions that may cause fluid film collapse and system failure. 
In steady-state operation, condensation may be an issue at low speeds even though 
GFBs usually heat up to a point where the saturated vapor pressure is never reached. 
To generally make GFBs less prone to vibrations at the cost of load-carrying capacity, 
the number and orientation of pads are identified as the most significant properties. 
Furthermore, the compliance and mobility of foil strips in combination with friction 
can be shown to play a decisive role in keeping the vibration amplitudes tolerable. 
Finally, the realistic addition of rotor unbalance does not invalidate previous findings 
but extends them to more complex scenarios with combined excitation mechanisms. 
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7 Conclusion 


7.2 


Perspective 


The findings of the present thesis give rise to a number of follow-up questions and 
thus provide interesting starting points for future research as briefly outlined below. 
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With regard to the fluid model, several simplifying assumptions are necessary 
in order to obtain a dynamical system and to limit the number of parameters. 
While cross-film gradients and phenomena such as wall-slip flow or turbulence 
most probably have no qualitative effect on the dynamics of the overall system, 
they should definitely be considered by rather quantitatively oriented studies. 


The regularized Coulomb friction fails in the prediction of sticking contacts, 
whereas more sophisticated models prove incompatible with standard methods. 
Consequently, future research that focuses on the foil structure should examine 
the applicability of numerical continuation for non-smooth dynamical systems. 


The considered Jeffcott-Laval rotor model makes no very specific assumptions 
and ignores the multi-disk design of most real-world refrigeration machinery. 
Also, it seems interesting to remove the restriction to symmetrical mode shapes 
and to investigate the amplitudes and frequencies of conical whirling motion. 


Since the lubrication gaps of GFBs are hardly accessible for measuring probes, 
the numerically predicted refrigerant behavior lacks experimental validation. 
Similarly, it might be interesting to compare the computed bifurcation diagrams 
and exemplary orbit shapes to data obtained from rotordynamic test benches. 


Finally, from a more general point of view, the elaborated modeling approach 
and solution strategy could be transferred not only to any kind of lubricant, 
but also to later GFB generations or alternative designs such as leaf-type GFBs. 


A Effect of Top Foil Deformations 


Top foil deformations are found to barely affect the rotordynamic behavior and can 
thus be neglected in studies focusing on qualitative tendencies of the overall system. 
However, it is quite interesting to look at the more and more wavy pressure fields as 
the top foil stiffness is lowered in Figs. A.1a and A.ıc compared to Fig. 6.3a for a 
three-pad GFB or in Figs. A.ıb and A.ıd compared to Fig. 6.5a for a single-pad GFB. 
For additional details about the underlying top foil model, which basically consists of 
segmented beams, the reader is referred to the Master’s thesis of the author [Lei15]. 
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(a) Pressure field (nominal top foil, Np = 3). (b) Pressure field (nominal top foil, Np = 1). 
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(c) Pressure field (softer top foil, Np = 3). (d) Pressure field (softer top foil, Np = 1). 


Figure A.1: Plots visualizing wavy fields during steady-state operation of three-pad or single-pad GFBs. 
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B Fluid Grid Independence Studies 


Prior to performing actual analyses, grid independence studies must be carried out 
in order to ensure the computational grid to be sufficiently fine for accurate results 
that do not change significantly upon further increases of the number of grid points. 
Taking the pressure field in a three-pad GFB from Fig. 6.3a with Nọ x Nz = 121 x 9, 
the plots in Figs. B.1a and B.1c show the found results on gradually coarser grids, 
which is done similarly for the single-pad GFB from Fig. 6.5a in Figs. B.1b and B.1d. 
Apparently, No x Nz = 25 x 5 already gives fairly good approximations of solutions. 
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(c) Pressure field (Np x Nz = 25 x 5, Np = 3). (d) Pressure field (Nọ x Nz = 25 x 5, Np = 1). 


Figure B.1: Plots visualizing approx. fields during steady-state operation of three-pad or single-pad GFBs. 
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Nomenclature 


Acronyms 

A/C Air Conditioning 

ACM Air Cycle Machine 

BVP Boundary Value Problem 

CFD Computational Fluid Dynamics 
DAE Differential-Algebraic Equation 
ECS Environmental Control System 
FDM Finite-Difference Method 

FEM Finite-Element Method 

FVM Finite-Volume Method 

GFB Gas Foil Bearing 

HB Hopf Bifurcation 

LP Limit Point 

LPC Limit Point of Cycles 

NSB Neimark-Sacker Bifurcation 
ODE Ordinary Differential Equation 
PDE Partial Differential Equation 

V Vapor (Region) 

VCM Vapor Cycle Machine 

V-L Vapor-Liquid (Region) 
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Nomenclature 


Greek Letters 


Symbol 
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Description Unit 
Compressibility coefficient field 1/K 
Normalized compressibility coefficient field — 
Combined heat transfer coefficient W/(m?-K) 
Rotor term of heat transfer coefficient Bo W/(m?:K) 


Structure term of heat transfer coefficient By W/(m? - K) 


Shaft journal attitude angle 
Nondimensional shaft journal attitude angle 
Rotor damping decay rate 

Nondimensional shaft journal eccentricity 
Nondimensional rotor mass eccentricity 
Small geometric parameter 

Vertical shaft journal position 

Vertical rotor disk position 

Vertical rotor disk mass position 
Nondimensional thermal management number 
Fluid temperature field 

Normalized fluid temperature field 

Fluid temperature (normal conditions) 
Characteristic fluid temperature difference 
Ambient fluid temperature 

Critical fluid temperature 

Nondimensional geometry number 


Nondimensional bearing number 


i-th bearing number on fixed-point solution branch 


A A A A 


Nomenclature 


Bearing number at Hopf bifurcation = (p. 71) 
Bearing number at limit point of cycles — (p. 74) 
Bearing number at Neimark-Sacker bifurcation — (p. 96) 
i-th bearing number on periodic solution branch = (p. 71) 
i-th frequency on periodic solution branch _ (p. 74) 
Volume viscosity field Pa:s (p. 27) 
Normalized volume viscosity field = (p. 28) 
n-th eigenvalue of JR = (p. 70) 
n-th eigenvalue of JR’ (Floquet multiplier) — (p. 73) 
Dynamic viscosity field Pa:s (p. 25) 
Normalized dynamic viscosity field = (p. 28) 
Dynamic viscosity (normal conditions) Pa-s (p. 26) 
Characteristic viscosity Pa:s (p. 28) 
Foil friction coefficient — (p. 52) 
Poisson’s ratio of foil material = (p. 47) 
Rotor mass distribution — (p. 16) 
Horizontal shaft journal position m (p. 13) 
Horizontal rotor disk position m (p. 17) 
Horizontal rotor disk mass position m (p. 17) 
Nondimensional fluid force (see F;,m) — (p. 57) 
Nondimensional horizontal bearing force = (p. 58) 
Nondimensional vertical bearing force = (p. 58) 
Nondimensional friction force (see fim) = (p. 57) 
Fluid density field kg/m? (p- 20) 
Normalized fluid density field = (p. 28) 
Fluid density (normal conditions) kg/m? (p. 56) 
Characteristic fluid density kg/m? (p- 28) 
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Nomenclature 


108 


Ambient fluid density 

Density of foil material 

Thickness of foil material 
Nondimensional time coordinate 

i-th period on periodic solution branch 
Nondimensional initial time 

Heat flux vector field 

Circumferential heat flux field 
Normalized circumferential heat flux field 
Radial heat flux field 

Normalized radial heat flux field 

Axial heat flux field 

Normalized axial heat flux field 

Angular coordinate 

Angular grid spacing 

Discretized angular coordinate 

Pad sector angle 

Angular position of n-th gap 

Bump sector angle 

Angular position of m-th bump on n-th pad 
Phase space of dynamical system 

General limit set in phase space Q 
Attracting limit set in phase space Q 
Repelling limit set in phase space Q 
Nondimensional bump frequency number 
Domain of n-th boundary value problem 


Boundary of n-th boundary value problem 


Nomenclature 


Mi Inflow region of boundary 0.2, — (p. 41) 
dOn, o Outflow region of boundary 00), — (p. 41) 
Mus Symmetry region of boundary 0), = (p. 41) 
OR Nondimensional rotor frequency number — (p. 58) 
w Angular velocity vector 1/s (p. 14) 
wo Angular frequency 1/s (p- 14) 
WB Bump natural angular frequency 1/s (p. 51) 
WR Rotor natural angular frequency 1/s (p. 18) 


Roman Letters 


Symbol Description Unit 

Ay Nondimensional compressibility coefficient field _ (p. 56) 
B Arbitrary fluid boundary point = (p. 14) 
Bg Nondimensional undeformed bump half-length — (p. 57) 
bg Undeformed bump half-length m (p. 46) 
BVP, Statement of n-th fluid BVP = (p. 44) 
C Nominal radial clearance m (p. 11) 
Co Characteristic specific heat capacity J/(kg -K) (p. 36) 
Cv Nondimensional specific heat capacity field — (p. 56) 
Cy Specific heat capacity field J/ (kg: K) (p. 23) 
Cv Normalized specific heat capacity field = (p. 36) 
C Specific heat capacity (normal conditions) J/(kg-K) (p. 26) 
D Nondimensional fluid density field — (p. 56) 
Dij Discretized fluid density — (p. 60) 
DR Nondimensional rotor damping number u (p. 58) 
AR Rotor damping coefficient N-s/m (p. 16) 
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Nomenclature 
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Specific internal energy field 

Shaft journal eccentricity 

Specific internal energy (normal conditions) 
Circumferential unit vector (rotor) 

Axial unit vector (rotor) 

Vertical unit vector (rotor) 

Horizontal unit vector (rotor) 

Young’s modulus of foil material 

Radial unit vector (rotor) 

Rotor mass eccentricity 

Circumferential unit vector (fluid/structure) 
Radial unit vector (fluid/structure) 

Axial unit vector (fluid/structure) 
Characteristic Eckert number 

Origin of coordinates (fluid /structure) 
Body force vector field 

Vertical bearing force 

Horizontal bearing force 

Fluid force for m-th bump on n-th pad 
Friction force for m-th bump on n-th pad 
Circumferential body force field 
Normalized circumferential body force field 
Radial body force field 

Normalized radial body force field 

Axial body force field 

Normalized axial body force field 


FDM right-hand side of Reynolds equation 


Nomenclature 


FDMr FDM right-hand side of temperature equation _ (p. 61) 
Fr Characteristic Froude number = (p. 29) 
G Nondimensional gravitational loading number = (p- 58) 
g Gravitational acceleration m/s* (p. 16) 
H Nondimensional fluid film thickness field — (p. 55) 
h Fluid film thickness field m (p. 14) 
hy Fluid film thickness field (rigid bearing) m (p. 14) 
ho Permissible fluid film thickness m (p. 15) 
Ag Nondimensional undeformed bump height — (p. 57) 
hg Undeformed bump height m (p. 46) 
H; Discretized fluid film thickness — (p. 59) 
Amin Minimum fluid film thickness m (p. 15) 
i Angular grid coordinate — (p. 59) 
j Axial grid coordinate = (p. 59) 
Jn Jacobian matrix for n-th pad = (p. 49) 
Jam Jacobian entry for m-th bump on n-th pad — (p. 49) 
je i-th Jacobian matrix on fixed-point solution branch — (p. 70) 
I: i-th Jacobian matrix on periodic solution branch ze (p. 73) 
kg Bump stiffness N/m (p. 47) 
kr Rotor shaft stiffness N/m (p. 16) 
Kz Nondimensional bristle stiffness number — (p. 58) 
kz Bristle stiffness N/m (p. 53) 
Kn Characteristic Knudsen number — (p. 20) 
L Axial bearing length m (p. 10) 
Llo Length scale of bearing dimensions m (p. 11) 
la Length scale of mean free path m (p. 20) 
En Lagrangian for n-th pad J (p. 50) 
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Nomenclature 


112 


Center of the shaft journal 

Center of the rotor disk 

Center of mass of the rotor disk 
Nondimensional bump mass number 
Bump mass 

Nondimensional rotor mass number 
Rotor mass 

i-th monodromy matrix on periodic solution branch 
Problem dimension 

Rotational speed 

Number of circumferential grid points 
Number of bumps per pad 

Number of internal nodes 


Number of pads 


Number of axial grid points 

Origin of coordinates (rotor) 

cos-like autonomization oscillator 
sin-like autonomization oscillator 
Nondimensional fluid pressure field 
Fluid pressure field 

Normalized fluid pressure field 

Fluid pressure (normal conditions) 
Nondimensional ambient fluid pressure 
Ambient fluid pressure 

Discretized fluid pressure 
Nondimensional foil structure deformation field 


Foil structure deformation field 


Nomenclature 


Discretized foil structure deformation — 
Generalized force for m-th bump on n-th pad N 
Shaft journal radius m 


Vector field of dynamical system = 


Shaft journal position vector m 
Rotor disk position vector m 
Rotor disk mass position vector m 
Specific gas constant J/ (kg: K) 


Characteristic Reynolds number — 
State vector — 
i-th state vector on fixed-point solution branch = 
i-th state vector on periodic solution branch = 
Perturbed state vector = 
Initial state vector = 
Fluid state vector = 
Rotor state vector ur 
Structure state vector = 
Nondimensional fluid temperature field = 
Time coordinate S 
Normalized time coordinate = 
Characteristic time scale S 
Nondimensional ambient fluid temperature = 


Discretized fluid temperature = 


Flow velocity vector field m/s 
Friction regularization velocity m/s 
Circumferential flow velocity field m/s 


Normalized circumferential flow velocity field — 
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Nomenclature 


114 


Circumferential bulk flow velocity field 

Couette term of flow velocity field 7, 

Poiseuille term of flow velocity field üx 

Radial flow velocity field 

Normalized radial flow velocity field 

Axial flow velocity field 

Normalized axial flow velocity field 

Axial bulk flow velocity field 

Nondimensional dynamic viscosity field 
Discretized dynamic viscosity 

Nondimensional generalized vapor quality field 
Generalized vapor quality field 

Generalized vapor quality (normal conditions) 
Nondimensional horizontal shaft journal position 
Circumferential coordinate (fluid/structure) 
Normalized circumferential coordinate 
Nondimensional generalized coordinate (see X, m) 
Generalized coordinate of m-th bump on n-th pad 
Nondimensional horizontal rotor disk position 
Nondimensional vertical shaft journal position 
Radial coordinate (fluid/structure) 


Normalized radial coordinate 


Nondimensional apex displacement (see yn,m) 
Apex displacement of m-th bump on n-th pad 
Nondimensional vertical rotor disk position 
Nondimensional axial coordinate 


Axial coordinate (fluid/structure) 


Nomenclature 


Axial grid spacing 

Normalized axial coordinate 

Discretized axial coordinate 
Nondimensional bristle deflection (see Zn,m) 


Bristle deflection of m-th bump on n-th pad 
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